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

      1: #ifndef included_ALE_LabelSifter_hh
      2: #define included_ALE_LabelSifter_hh

      4: #include <iostream>

      6: #ifndef  included_ALE_hh
      7: #include <sieve/ALE.hh>
      8: #endif

     10: namespace ALE {
     11:   namespace NewSifterDef {
     12:     // Defines the traits of a sequence representing a subset of a multi_index container Index_.
     13:     // A sequence defines output (input in std terminology) iterators for traversing an Index_ object.
     14:     // Upon dereferencing values are extracted from each result record using a ValueExtractor_ object.
     15:     template <typename Index_, typename ValueExtractor_>
     16:     struct IndexSequenceTraits {
     17:       typedef Index_ index_type;
     18:       class iterator_base {
     19:       public:
     20:         // Standard iterator typedefs
     21:         typedef ValueExtractor_                        extractor_type;
     22:         typedef std::input_iterator_tag                iterator_category;
     23:         typedef typename extractor_type::result_type   value_type;
     24:         typedef int                                    difference_type;
     25:         typedef value_type*                            pointer;
     26:         typedef value_type&                            reference;

     28:         // Underlying iterator type
     29:         typedef typename index_type::iterator          itor_type;
     30:       protected:
     31:         // Underlying iterator
     32:         itor_type      _itor;
     33:         // Member extractor
     34:         extractor_type _ex;
     35:       public:
     36:         iterator_base(itor_type itor) {
     37:           this->_itor = itor_type(itor);
     38:         };
     39:         virtual ~iterator_base() {};
     40:         virtual bool              operator==(const iterator_base& iter) const {return this->_itor == iter._itor;};
     41:         virtual bool              operator!=(const iterator_base& iter) const {return this->_itor != iter._itor;};
     42:         // FIX: operator*() should return a const reference, but it won't compile that way, because _ex() returns const value_type
     43:         virtual const value_type  operator*() const {return _ex(*(this->_itor));};
     44:       };// class iterator_base
     45:       class iterator : public iterator_base {
     46:       public:
     47:         // Standard iterator typedefs
     48:         typedef typename iterator_base::iterator_category  iterator_category;
     49:         typedef typename iterator_base::value_type         value_type;
     50:         typedef typename iterator_base::extractor_type     extractor_type;
     51:         typedef typename iterator_base::difference_type    difference_type;
     52:         typedef typename iterator_base::pointer            pointer;
     53:         typedef typename iterator_base::reference          reference;
     54:         // Underlying iterator type
     55:         typedef typename iterator_base::itor_type          itor_type;
     56:       public:
     57:         iterator(const itor_type& itor) : iterator_base(itor) {};
     58:         virtual ~iterator() {};
     59:         //
     60:         virtual iterator   operator++() {++this->_itor; return *this;};
     61:         virtual iterator   operator++(int n) {iterator tmp(this->_itor); ++this->_itor; return tmp;};
     62:       };// class iterator
     63:     }; // struct IndexSequenceTraits

     65:     template <typename Index_, typename ValueExtractor_>
     66:     struct ReversibleIndexSequenceTraits {
     67:       typedef IndexSequenceTraits<Index_, ValueExtractor_> base_traits;
     68:       typedef typename base_traits::iterator_base   iterator_base;
     69:       typedef typename base_traits::iterator        iterator;
     70:       typedef typename base_traits::index_type      index_type;

     72:       // reverse_iterator is the reverse of iterator
     73:       class reverse_iterator : public iterator_base {
     74:       public:
     75:         // Standard iterator typedefs
     76:         typedef typename iterator_base::iterator_category  iterator_category;
     77:         typedef typename iterator_base::value_type         value_type;
     78:         typedef typename iterator_base::extractor_type     extractor_type;
     79:         typedef typename iterator_base::difference_type    difference_type;
     80:         typedef typename iterator_base::pointer            pointer;
     81:         typedef typename iterator_base::reference          reference;
     82:         // Underlying iterator type
     83:         typedef typename iterator_base::itor_type          itor_type;
     84:       public:
     85:         reverse_iterator(const itor_type& itor) : iterator_base(itor) {};
     86:         virtual ~reverse_iterator() {};
     87:         //
     88:         virtual reverse_iterator     operator++() {--this->_itor; return *this;};
     89:         virtual reverse_iterator     operator++(int n) {reverse_iterator tmp(this->_itor); --this->_itor; return tmp;};
     90:       };
     91:     }; // class ReversibleIndexSequenceTraits

     93:     //
     94:     // Arrow & ArrowContainer definitions
     95:     //
     96:     template<typename Source_, typename Target_>
     97:     struct  Arrow { //: public ALE::def::Arrow<Source_, Target_, Color_> {
     98:       typedef Arrow   arrow_type;
     99:       typedef Source_ source_type;
    100:       typedef Target_ target_type;
    101:       source_type source;
    102:       target_type target;
    103:       Arrow(const source_type& s, const target_type& t) : source(s), target(t) {};
    104:       // Flipping
    105:       template <typename OtherSource_, typename OtherTarget_>
    106:       struct rebind {
    107:         typedef Arrow<OtherSource_, OtherTarget_> type;
    108:       };
    109:       struct flip {
    110:         typedef Arrow<target_type, source_type> type;
    111:         type arrow(const arrow_type& a) { return type(a.target, a.source);};
    112:       };

    114:       // Printing
    115:       friend std::ostream& operator<<(std::ostream& os, const Arrow& a) {
    116:         os << a.source << " ----> " << a.target;
    117:         return os;
    118:       }

    120:       // Arrow modifiers
    121:       struct sourceChanger {
    122:         sourceChanger(const source_type& newSource) : _newSource(newSource) {};
    123:         void operator()(arrow_type& a) {a.source = this->_newSource;}
    124:       private:
    125:         source_type _newSource;
    126:       };

    128:       struct targetChanger {
    129:         targetChanger(const target_type& newTarget) : _newTarget(newTarget) {};
    130:         void operator()(arrow_type& a) { a.target = this->_newTarget;}
    131:       private:
    132:         const target_type _newTarget;
    133:       };
    134:     };// struct Arrow

    137:     template<typename Source_, typename Target_>
    138:     struct ArrowContainerTraits {
    139:     public:
    140:       //
    141:       // Encapsulated types
    142:       //
    143:       typedef Arrow<Source_,Target_>           arrow_type;
    144:       typedef typename arrow_type::source_type source_type;
    145:       typedef typename arrow_type::target_type target_type;
    146:       // Index tags
    147:       struct                                   sourceTargetTag{};
    148:       struct                                   targetSourceTag{};

    150:       // Sequence traits and sequence types
    151:       template <typename Index_, typename Key_, typename SubKey_, typename ValueExtractor_>
    152:       class ArrowSequence {
    153:         // ArrowSequence implements ReversibleIndexSequencTraits with Index_ and ValueExtractor_ types.
    154:         // A Key_ object and an optional SubKey_ object are used to extract the index subset.
    155:       public:
    156:         typedef ReversibleIndexSequenceTraits<Index_, ValueExtractor_>  traits;
    157:         //typedef source_type                                             source_type;
    158:         //typedef target_type                                             target_type;
    159:         //typedef arrow_type                                              arrow_type;
    160:         //
    161:         typedef Key_                                                    key_type;
    162:         typedef SubKey_                                                 subkey_type;
    163:       protected:
    164:         typename traits::index_type&                                    _index;
    165:         key_type                                                  key;
    166:         subkey_type                                               subkey;
    167:         bool                                                      useSubkey;
    168:       public:
    169:         // Need to extend the inherited iterators to be able to extract arrow color
    170:         class iterator : public traits::iterator {
    171:         public:
    172:           iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
    173:           virtual const source_type& source() const {return this->_itor->source;};
    174:           virtual const target_type& target() const {return this->_itor->target;};
    175:           virtual const arrow_type&  arrow()  const {return *(this->_itor);};
    176:         };
    177:         class reverse_iterator : public traits::reverse_iterator {
    178:         public:
    179:           reverse_iterator(const typename traits::reverse_iterator::itor_type& itor) : traits::reverse_iterator(itor) {};
    180:           virtual const source_type& source() const {return this->_itor->source;};
    181:           virtual const target_type& target() const {return this->_itor->target;};
    182:           virtual const arrow_type&  arrow()  const {return *(this->_itor);};
    183:         };
    184:       public:
    185:         //
    186:         // Basic ArrowSequence interface
    187:         //
    188:         ArrowSequence(const ArrowSequence& seq) : _index(seq._index), key(seq.key), subkey(seq.subkey), useSubkey(seq.useSubkey) {};
    189:         ArrowSequence(typename traits::index_type& index, const key_type& k) :
    190:           _index(index), key(k), subkey(subkey_type()), useSubkey(0) {};
    191:         ArrowSequence(typename traits::index_type& index, const key_type& k, const subkey_type& kk) :
    192:           _index(index), key(k), subkey(kk), useSubkey(1){};
    193:         virtual ~ArrowSequence() {};

    195:         void setKey(const key_type& key) {this->key = key;};
    196:         void setSubkey(const subkey_type& subkey) {this->subkey = subkey;};
    197:         void setUseSubkey(const bool& useSubkey) {this->useSubkey = useSubkey;};

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

    201:         virtual typename traits::index_type::size_type  size()  {
    202:           if (this->useSubkey) {
    203:             return this->_index.count(::boost::make_tuple(this->key,this->subkey));
    204:           } else {
    205:             return this->_index.count(::boost::make_tuple(this->key));
    206:           }
    207:         };

    209:         virtual iterator begin() {
    210:           if (this->useSubkey) {
    211:             return iterator(this->_index.lower_bound(::boost::make_tuple(this->key,this->subkey)));
    212:           } else {
    213:             return iterator(this->_index.lower_bound(::boost::make_tuple(this->key)));
    214:           }
    215:         };

    217:         virtual iterator end() {
    218:           if (this->useSubkey) {
    219:             return iterator(this->_index.upper_bound(::boost::make_tuple(this->key,this->subkey)));
    220:           } else {
    221:             return iterator(this->_index.upper_bound(::boost::make_tuple(this->key)));
    222:           }
    223:         };

    225:         virtual reverse_iterator rbegin() {
    226:           if (this->useSubkey) {
    227:             return reverse_iterator(--this->_index.upper_bound(::boost::make_tuple(this->key,this->subkey)));
    228:           } else {
    229:             return reverse_iterator(--this->_index.upper_bound(::boost::make_tuple(this->key)));
    230:           }
    231:         };

    233:         virtual reverse_iterator rend() {
    234:           if (this->useSubkey) {
    235:             return reverse_iterator(--this->_index.lower_bound(::boost::make_tuple(this->key,this->subkey)));
    236:           } else {
    237:             return reverse_iterator(--this->_index.lower_bound(::boost::make_tuple(this->key)));
    238:           }
    239:         };

    241:         template<typename ostream_type>
    242:         void view(ostream_type& os, const char* label = NULL){
    243:           if(label != NULL) {
    244:             os << "Viewing " << label << " sequence:" << std::endl;
    245:           }
    246:           os << "[";
    247:           for(iterator i = this->begin(); i != this->end(); i++) {
    248:             os << " (" << *i << ")";
    249:           }
    250:           os << " ]" << std::endl;
    251:         }
    252:       };// class ArrowSequence
    253:     };// class ArrowContainerTraits

    256:     // The specialized ArrowContainer types distinguish the cases of unique and multiple colors of arrows on
    257:     // for each (source,target) pair (i.e., a single arrow, or multiple arrows between each pair of points).
    258:     template<typename Source_, typename Target_, typename Alloc_ = ALE_ALLOCATOR<typename ArrowContainerTraits<Source_, Target_>::arrow_type> >
    259:     struct ArrowContainer {
    260:       // Define container's encapsulated types
    261:       typedef ArrowContainerTraits<Source_, Target_> traits;
    262:       // need to def arrow_type locally, since BOOST_MULTI_INDEX_MEMBER barfs when first template parameter starts with 'typename'
    263:       typedef typename traits::arrow_type                                   arrow_type;
    264:       typedef Alloc_ alloc_type;

    266:       // multi-index set type -- arrow set
    267:       typedef ::boost::multi_index::multi_index_container<
    268:         typename traits::arrow_type,
    269:         ::boost::multi_index::indexed_by<
    270:           ::boost::multi_index::ordered_unique<
    271:             ::boost::multi_index::tag<typename traits::sourceTargetTag>,
    272:             ::boost::multi_index::composite_key<
    273:               typename traits::arrow_type,
    274:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source),
    275:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target)
    276:             >
    277:           >,
    278:           ::boost::multi_index::ordered_unique<
    279:             ::boost::multi_index::tag<typename traits::targetSourceTag>,
    280:             ::boost::multi_index::composite_key<
    281:               typename traits::arrow_type,
    282:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target),
    283:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source)
    284:             >
    285:           >
    286:         >,
    287:         Alloc_
    288:       > set_type;
    289:       // multi-index set of arrow records
    290:       set_type set;

    292:       ArrowContainer() {};
    293:       ArrowContainer(const alloc_type& allocator) {this->set = set_type(typename set_type::ctor_args_list(), allocator);};
    294:     }; // class ArrowContainer
    295:   }; // namespace NewSifterDef

    297:   template<typename Source_, typename Target_, typename Alloc_ = ALE_ALLOCATOR<typename NewSifterDef::ArrowContainer<Source_, Target_>::traits::arrow_type> >
    298:   class LabelSifter { // class Sifter
    299:   public:
    300:     typedef struct {
    301:       typedef LabelSifter<Source_, Target_, Alloc_> graph_type;
    302:       // Encapsulated container types
    303:       typedef NewSifterDef::ArrowContainer<Source_, Target_, Alloc_>                 arrow_container_type;
    304:       // Types associated with records held in containers
    305:       typedef typename arrow_container_type::traits::arrow_type                      arrow_type;
    306:       typedef typename arrow_container_type::traits::source_type                     source_type;
    307:       typedef typename arrow_container_type::traits::target_type                     target_type;
    308:       // Convenient tag names
    309:       typedef typename arrow_container_type::traits::sourceTargetTag                 supportInd;
    310:       typedef typename arrow_container_type::traits::targetSourceTag                 coneInd;
    311:       typedef typename arrow_container_type::traits::sourceTargetTag                 arrowInd;
    312:       //
    313:       // Return types
    314:       //
    315:       typedef typename
    316:       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, source_type, source)>
    317:       arrowSequence;

    319:       // FIX: This is a temp fix to include addArrow into the interface; should probably be pushed up to ArrowSequence
    320:       struct coneSequence : public arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,coneInd>::type, target_type, source_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, source_type, source)> {
    321:       protected:
    322:         graph_type& _graph;
    323:       public:
    324:         typedef typename
    325:           arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,coneInd>::type, target_type, source_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, source_type, source)> base_type;
    326:         // Encapsulated types
    327:         typedef typename base_type::traits traits;
    328:         typedef typename base_type::iterator iterator;
    329:         typedef typename base_type::reverse_iterator reverse_iterator;
    330:         // Basic interface
    331:         coneSequence(const coneSequence& seq) : base_type(seq), _graph(seq._graph) {};
    332:           coneSequence(graph_type& graph, typename traits::index_type& index, const typename base_type::key_type& k) : base_type(index, k), _graph(graph){};
    333:             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) {};
    334:               virtual ~coneSequence() {};

    336:         // Extended interface
    337:         void addArrow(const arrow_type& a) {
    338:           // if(a.target != this->key) {
    339:           //               throw ALE::Exception("Arrow target mismatch in a coneSequence");
    340:           //             }
    341:           this->_graph.addArrow(a);
    342:         };
    343:         void addArrow(const source_type& s){
    344:           this->_graph.addArrow(arrow_type(s,this->key));
    345:         };

    347:         virtual bool contains(const source_type& s) {
    348:           // Check whether a given point is in the index
    349:           typename ::boost::multi_index::index<typename LabelSifter::traits::arrow_container_type::set_type,typename LabelSifter::traits::arrowInd>::type& index = ::boost::multi_index::get<typename LabelSifter::traits::arrowInd>(this->_graph._arrows.set);
    350:           return (index.find(::boost::make_tuple(s,this->key)) != index.end());
    351:         };
    352:       };// struct coneSequence

    354:       // FIX: This is a temp fix to include addArrow into the interface; should probably be pushed up to ArrowSequence
    355:       struct supportSequence : public arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,supportInd>::type, source_type, target_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, target_type, target)> {
    356:       protected:
    357:         graph_type& _graph;
    358:       public:
    359:         typedef typename
    360:           arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,supportInd>::type, source_type, target_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, target_type, target)> base_type;
    361:         // Encapsulated types
    362:         typedef typename base_type::traits traits;
    363:         typedef typename base_type::iterator iterator;
    364:         typedef typename base_type::iterator const_iterator;
    365:         typedef typename base_type::reverse_iterator reverse_iterator;
    366:         // Basic interface
    367:         supportSequence(const supportSequence& seq) : base_type(seq), _graph(seq._graph) {};
    368:         supportSequence(graph_type& graph, typename traits::index_type& index, const typename base_type::key_type& k) : base_type(index, k), _graph(graph){};
    369:         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) {};
    370:         virtual ~supportSequence() {};

    372:         // FIX: WARNING: (or a HACK?): we flip the arrow on addition here.
    373:         // Fancy interface
    374:         void addArrow(const typename arrow_type::flip::type& af) {
    375:           this->_graph.addArrow(af.target, af.source);
    376:         };
    377:         void addArrow(const target_type& t){
    378:           this->_graph.addArrow(arrow_type(this->key,t));
    379:         };
    380:       };// struct supportSequence

    382:       typedef std::set<source_type, std::less<source_type>, typename Alloc_::template rebind<source_type>::other> coneSet;
    383:       typedef ALE::array<source_type> coneArray;
    384:       typedef std::set<target_type, std::less<target_type>, typename Alloc_::template rebind<source_type>::other> supportSet;
    385:       typedef ALE::array<target_type> supportArray;
    386:     } traits;

    388:     template <typename OtherSource_, typename OtherTarget_>
    389:     struct rebind {
    390:       typedef LabelSifter<OtherSource_, OtherTarget_> type;
    391:     };

    393:     typedef Alloc_                           alloc_type;
    394:     typedef typename traits::source_type     source_type;
    395:     typedef typename traits::target_type     target_type;
    396:     typedef typename traits::coneSequence    coneSequence;
    397:     typedef typename traits::supportSequence supportSequence;
    398:     typedef std::set<int>                    capSequence;
    399:   public:
    400:     // Debug level
    401:     int _debug;
    402:     //protected:
    403:     typename traits::arrow_container_type _arrows;
    404:   protected:
    405:     MPI_Comm    _comm;
    406:     int         _commRank;
    407:     int         _commSize;
    408:     void __init(MPI_Comm comm) {
    409:       static PetscClassId sifterType = -1;
    410:       //const char        *id_name = ALE::getClassName<T>();
    411:       const char        *id_name = "LabelSifter";
    412:       PetscErrorCode     ierr;

    414:       if (sifterType < 0) {
    415:         PetscClassIdRegister(id_name,&sifterType);CHKERROR(ierr, "Error in MPI_Comm_rank");
    416:       }
    417:       this->_comm = comm;
    418:       MPI_Comm_rank(this->_comm, &this->_commRank);CHKERROR(ierr, "Error in MPI_Comm_rank");
    419:       MPI_Comm_size(this->_comm, &this->_commSize);CHKERROR(ierr, "Error in MPI_Comm_rank");
    420:       //ALE::restoreClassName<T>(id_name);
    421:     };
    422:     // We store these sequence objects to avoid creating them each query
    423:     Obj<typename traits::coneSequence> _coneSeq;
    424:     Obj<typename traits::supportSequence> _supportSeq;
    425:   public:
    426:     //
    427:     // Basic interface
    428:     //
    429:     LabelSifter(MPI_Comm comm = PETSC_COMM_SELF, const int& debug = 0) : _debug(debug) {
    430:       __init(comm);
    431:       this->_coneSeq    = new typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), typename traits::target_type());
    432:       this->_supportSeq = new typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), typename traits::source_type());
    433:     };
    434:     LabelSifter(MPI_Comm comm, Alloc_& allocator, const int& debug) : _debug(debug), _arrows(allocator) {
    435:       __init(comm);
    436:       this->_coneSeq    = new typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), typename traits::target_type());
    437:       this->_supportSeq = new typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), typename traits::source_type());
    438:     };
    439:     virtual ~LabelSifter() {};
    440:     //
    441:     // Query methods
    442:     //
    443:     int         debug()    const {return this->_debug;};
    444:     void        setDebug(const int debug) {this->_debug = debug;};
    445:     MPI_Comm    comm()     const {return this->_comm;};
    446:     int         commSize() const {return this->_commSize;};
    447:     int         commRank() const {return this->_commRank;}

    449:     // FIX: should probably have cone and const_cone etc, since arrows can be modified through an iterator (modifyColor).
    450:     Obj<typename traits::arrowSequence>
    451:     arrows(const typename traits::source_type& s, const typename traits::target_type& t) {
    452:       return typename traits::arrowSequence(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s, t);
    453:     };
    454:     Obj<typename traits::arrowSequence>
    455:     arrows(const typename traits::source_type& s) {
    456:       return typename traits::arrowSequence(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s);
    457:     };
    458: #ifdef SLOW
    459:     Obj<typename traits::coneSequence>
    460:     cone(const typename traits::target_type& p) {
    461:       return typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);
    462:     };
    463: #else
    464:     const Obj<typename traits::coneSequence>&
    465:     cone(const typename traits::target_type& p) {
    466:       this->_coneSeq->setKey(p);
    467:       this->_coneSeq->setUseSubkey(false);
    468:       return this->_coneSeq;
    469:     };
    470: #endif
    471:     template<class InputSequence>
    472:     Obj<typename traits::coneSet>
    473:     cone(const Obj<InputSequence>& points) {
    474:       Obj<typename traits::coneSet> cone = typename traits::coneSet();

    476:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
    477:         const Obj<typename traits::coneSequence>& pCone = this->cone(*p_itor);
    478:         cone->insert(pCone->begin(), pCone->end());
    479:       }
    480:       return cone;
    481:     }
    482:     int getConeSize(const typename traits::target_type& p) {
    483:       return this->cone(p)->size();
    484:     };
    485:     template<typename PointCheck>
    486:     bool coneContains(const typename traits::target_type& p, const PointCheck& checker) {
    487:       typename traits::coneSequence cone(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);

    489:       for(typename traits::coneSequence::iterator c_iter = cone.begin(); c_iter != cone.end(); ++c_iter) {
    490:         if (checker(*c_iter, p)) return true;
    491:       }
    492:       return false;
    493:     }
    494:     template<typename PointProcess>
    495:     void coneApply(const typename traits::target_type& p, PointProcess& processor) {
    496:       typename traits::coneSequence cone(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);

    498:       for(typename traits::coneSequence::iterator c_iter = cone.begin(); c_iter != cone.end(); ++c_iter) {
    499:         processor(*c_iter, p);
    500:       }
    501:     }
    502: #ifdef SLOW
    503:     Obj<typename traits::supportSequence>
    504:     support(const typename traits::source_type& p) {
    505:       return typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);
    506:     };
    507: #else
    508:     const Obj<typename traits::supportSequence>&
    509:     support(const typename traits::source_type& p) {
    510:       this->_supportSeq->setKey(p);
    511:       this->_supportSeq->setUseSubkey(false);
    512:       return this->_supportSeq;
    513:     };
    514: #endif
    515:     template<class InputSequence>
    516:     Obj<typename traits::supportSet>
    517:     support(const Obj<InputSequence>& points){
    518:       Obj<typename traits::supportSet> supp = typename traits::supportSet();
    519:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
    520:         const Obj<typename traits::supportSequence>& pSupport = this->support(*p_itor);
    521:         supp->insert(pSupport->begin(), pSupport->end());
    522:       }
    523:       return supp;
    524:     }
    525:     template<typename PointCheck>
    526:     bool supportContains(const typename traits::source_type& p, const PointCheck& checker) {
    527:       typename traits::supportSequence support(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);

    529:       for(typename traits::supportSequence::iterator s_iter = support.begin(); s_iter != support.end(); ++s_iter) {
    530:         if (checker(*s_iter, p)) return true;
    531:       }
    532:       return false;
    533:     }
    534:     template<typename PointProcess>
    535:     void supportApply(const typename traits::source_type& p, PointProcess& processor) {
    536:       typename traits::supportSequence support(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);

    538:       for(typename traits::supportSequence::iterator s_iter = support.begin(); s_iter != support.end(); ++s_iter) {
    539:         processor(*s_iter, p);
    540:       }
    541:     }

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

    547:       if(label != NULL) {
    548:         os << "["<<rank<<"]Viewing LabelSifter '" << label << "':" << std::endl;
    549:       }
    550:       else {
    551:         os << "["<<rank<<"]Viewing a LabelSifter:" << std::endl;
    552:       }
    553:       os << "'raw' arrow set:" << std::endl;
    554:       for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++) {
    555:         os << *ai << std::endl;
    556:       }
    557:     }
    558:     // A parallel viewer
    561:     PetscErrorCode view(const char* label = NULL, bool raw = false){
    563:       ostringstream txt;
    565:       if(this->_debug) {
    566:         std::cout << "viewing a LabelSifter, comm = " << this->comm() << ", PETSC_COMM_SELF = " << PETSC_COMM_SELF << ", commRank = " << this->commRank() << std::endl;
    567:       }
    568:       if(label != NULL) {
    569:         PetscPrintf(this->comm(), "viewing LabelSifter: '%s'\n", label);
    570:       } else {
    571:         PetscPrintf(this->comm(), "viewing a LabelSifter: \n");
    572:       }
    573:       if(!raw) {
    574:         ostringstream txt;
    575:         if(this->commRank() == 0) {
    576:           txt << "cap --> base:\n";
    577:         }
    578:         if(_arrows.set.empty()) {
    579:           txt << "[" << this->commRank() << "]: empty" << std::endl;
    580:         }
    581:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++) {
    582:           txt << "[" << this->commRank() << "]: " << ai->source << "---->" << ai->target << std::endl;
    583:         }
    584:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    585:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    586:       }
    587:       else { // if(raw)
    588:         ostringstream txt;
    589:         if(this->commRank() == 0) {
    590:           txt << "'raw' arrow set:" << std::endl;
    591:         }
    592:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++)
    593:         {
    594:           typename traits::arrow_type arr = *ai;
    595:           txt << "[" << this->commRank() << "]: " << arr << std::endl;
    596:         }
    597:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    598:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    599:       }// if(raw)

    601:       return(0);
    602:     }
    603:   public:
    604:     //
    605:     // Lattice queries
    606:     //
    607:     template<class targetInputSequence>
    608:     Obj<typename traits::coneSequence> meet(const Obj<targetInputSequence>& targets);
    609:     // unimplemented
    610:     template<class sourceInputSequence>
    611:     Obj<typename traits::coneSequence> join(const Obj<sourceInputSequence>& sources);
    612:   public:
    613:     //
    614:     // Structural manipulation
    615:     //
    616:     void clear() {
    617:       this->_arrows.set.clear();
    618:     };
    619:     // This is necessary to work with Completion right now
    620:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q, const int dummy) {
    621:       this->addArrow(p, q);
    622:     };
    623:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q) {
    624:       this->addArrow(typename traits::arrow_type(p, q));
    625:       //std::cout << "Added " << arrow_type(p, q);
    626:     };
    627:     virtual void addArrow(const typename traits::arrow_type& a) {
    628:       this->_arrows.set.insert(a);
    629:     };
    630:     virtual void removeArrow(const typename traits::arrow_type& a) {
    631:       // First, produce an arrow sequence for the given source, target combination.
    632:       typename traits::arrowSequence::traits::index_type& arrowIndex =
    633:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
    634:       typename traits::arrowSequence::traits::index_type::iterator i,ii,j;
    635:       i = arrowIndex.lower_bound(::boost::make_tuple(a.source,a.target));
    636:       ii = arrowIndex.upper_bound(::boost::make_tuple(a.source, a.target));
    637:       if (this->_debug) {
    638:         std::cout << "removeArrow: attempting to remove arrow:" << a << std::endl;
    639:         std::cout << "removeArrow: candidate arrows are:" << std::endl;
    640:       }
    641:       for(j = i; j != ii; j++) {
    642:         if (this->_debug) {
    643:           std::cout << " " << *j;
    644:         }
    645:         // Find the arrow of right color and remove it
    646:         if (this->_debug) {
    647:           std::cout << std::endl << "removeArrow: found:" << *j << std::endl;
    648:         }
    649:         arrowIndex.erase(j);
    650:         break;
    651:       }
    652:     };

    654:     void addCone(const typename traits::source_type& source, const typename traits::target_type& target){
    655:       this->addArrow(source, target);
    656:     };
    657:     template<class sourceInputSequence>
    658:     void
    659:     addCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target){
    660:       if (this->_debug > 1) {std::cout << "Adding a cone " << std::endl;}
    661:       for(typename sourceInputSequence::iterator iter = sources->begin(); iter != sources->end(); ++iter) {
    662:         if (this->_debug > 1) {std::cout << "Adding arrow from " << *iter << " to " << target << std::endl;}
    663:         this->addArrow(*iter, target);
    664:       }
    665:     }
    666:     void clearCone(const typename traits::target_type& t) {
    667:       // Use the cone sequence types to clear the cone
    668:       typename traits::coneSequence::traits::index_type& coneIndex =
    669:         ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set);
    670:       typename traits::coneSequence::traits::index_type::iterator i, ii, j;
    671:       if (this->_debug > 20) {
    672:         std::cout << "clearCone: removing cone over " << t;
    673:         std::cout << std::endl;
    674:         const Obj<typename traits::coneSequence>& cone = this->cone(t);
    675:         std::cout << "[";
    676:         for(typename traits::coneSequence::iterator ci = cone->begin(); ci != cone->end(); ci++) {
    677:           std::cout << "  " << ci.arrow();
    678:         }
    679:         std::cout << "]" << std::endl;
    680:       }
    681:       i = coneIndex.lower_bound(::boost::make_tuple(t));
    682:       ii = coneIndex.upper_bound(::boost::make_tuple(t));
    683:       coneIndex.erase(i,ii);
    684:     }// clearCone()

    686:     void clearSupport(const typename traits::source_type& s) {
    687:       // Use the cone sequence types to clear the cone
    688:       typename
    689:         traits::supportSequence::traits::index_type& suppIndex = ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set);
    690:       typename traits::supportSequence::traits::index_type::iterator i, ii, j;
    691:       i = suppIndex.lower_bound(::boost::make_tuple(s));
    692:       ii = suppIndex.upper_bound(::boost::make_tuple(s));
    693:       suppIndex.erase(i,ii);
    694:     }
    695:     void setCone(const typename traits::source_type& source, const typename traits::target_type& target){
    696:       this->clearCone(target); this->addCone(source, target);
    697:     }
    698:     template<class sourceInputSequence>
    699:     void setCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target) {
    700:       this->clearCone(target); this->addCone(sources, target);
    701:     }
    702:     template<class targetInputSequence>
    703:     void addSupport(const typename traits::source_type& source, const Obj<targetInputSequence >& targets) {
    704:       if (this->_debug > 1) {std::cout << "Adding a support " << std::endl;}
    705:       for(typename targetInputSequence::iterator iter = targets->begin(); iter != targets->end(); ++iter) {
    706:         if (this->_debug > 1) {std::cout << "Adding arrow from " << source << " to " << *iter << std::endl;}
    707:         this->addArrow(source, *iter);
    708:       }
    709:     }
    710:     template<typename Sifter_, typename AnotherSifter_>
    711:     void add(const Obj<Sifter_>& cbg, const Obj<AnotherSifter_>& baseRestriction = NULL) {
    712:       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);
    713:       bool baseRestrict = !baseRestriction.isNull();

    715:       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) {
    716:         if (baseRestrict) {
    717:           if (!baseRestriction->getSupportSize(a_iter->target) && !baseRestriction->getConeSize(a_iter->target)) continue;
    718:         }
    719:         this->addArrow(*a_iter);
    720:       }
    721:     }
    722:     template<typename Sifter_, typename AnotherSifter_, typename Renumbering_>
    723:     void add(const Obj<Sifter_>& cbg, const Obj<AnotherSifter_>& baseRestriction, Renumbering_& renumbering) {
    724:       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);

    726:       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) {
    727:         if (renumbering.find(a_iter->target) == renumbering.end()) continue;
    728:         target_type target = renumbering[a_iter->target];

    730:         if (!baseRestriction->getSupportSize(target) && !baseRestriction->getConeSize(target)) continue;
    731:         this->addArrow(a_iter->source, target);
    732:       }
    733:     }
    734:     template<typename Labeling, typename AnotherSifter>
    735:     void relabel(Labeling& relabeling, AnotherSifter& newLabel) {
    736:       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);

    738:       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) {
    739:         const typename traits::target_type newTarget = relabeling.restrictPoint(a_iter->target)[0];

    741:         newLabel.addArrow(a_iter->source, newTarget);
    742:       }
    743:     }

    745:     int size() const {return _arrows.set.size();};
    746:     int getCapSize() const {
    747:       std::set<source_type> cap;
    748:       for(typename traits::arrow_container_type::set_type::iterator a_iter = _arrows.set.begin(); a_iter != _arrows.set.end(); ++a_iter) {
    749:         cap.insert(a_iter->source);
    750:       }
    751:       return cap.size();
    752:     };
    753:     capSequence cap() const {
    754:       std::set<source_type> cap;
    755:       for(typename traits::arrow_container_type::set_type::iterator a_iter = _arrows.set.begin(); a_iter != _arrows.set.end(); ++a_iter) {
    756:         cap.insert(a_iter->source);
    757:       }
    758:       return cap;
    759:     };
    760:     int getBaseSize() const {
    761:       std::set<target_type> base;
    762:       for(typename traits::arrow_container_type::set_type::iterator a_iter = _arrows.set.begin(); a_iter != _arrows.set.end(); ++a_iter) {
    763:         base.insert(a_iter->target);
    764:       }
    765:       return base.size();
    766:     };
    767:   public: // Compatibility with fixed storage variants
    768:     typedef Interval<target_type> chart_type;
    769:     chart_type& getChart() {static chart_type chart(0, 0); return chart;}
    770:     template<typename chart_type>
    771:     void setChart(const chart_type& chart) {}
    772:     void setConeSize(target_type p, int s) {}
    773:     void setSupportSize(source_type p, int s) {}
    774:     void allocate() {}
    775:     void recalculateLabel() {}
    776:   }; // class LabelSifter

    778:   class LabelSifterSerializer {
    779:   public:
    780:     template<typename LabelSifter>
    781:     static void writeLabel(std::ofstream& fs, LabelSifter& label) {
    782:       if (label.commRank() == 0) {
    783:         // Write local
    784:         fs << label._arrows.set.size() << std::endl;
    785:         for(typename LabelSifter::traits::arrow_container_type::set_type::iterator ai = label._arrows.set.begin(); ai != label._arrows.set.end(); ai++) {
    786:           fs << ai->source << " " << ai->target << std::endl;
    787:         }
    788:         // Receive and write remote
    789:         for(int p = 1; p < label.commSize(); ++p) {
    790:           PetscInt       size;
    791:           PetscInt      *arrows;
    792:           MPI_Status     status;

    795:           MPI_Recv(&size, 1, MPIU_INT, p, 1, label.comm(), &status);CHKERRXX(ierr);
    796:           fs << size << std::endl;
    797:           PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    798:           MPI_Recv(arrows, size*2, MPIU_INT, p, 1, label.comm(), &status);CHKERRXX(ierr);
    799:           for(PetscInt a = 0; a < size; ++a) {
    800:             fs << arrows[a*2+0] << " " << arrows[a*2+1] << std::endl;
    801:           }
    802:           PetscFree(arrows);CHKERRXX(ierr);
    803:         }
    804:       } else {
    805:         // Send remote
    806:         PetscInt       size = label._arrows.set.size();
    807:         PetscInt       a    = 0;
    808:         PetscInt      *arrows;

    811:         MPI_Send(&size, 1, MPIU_INT, 0, 1, label.comm());CHKERRXX(ierr);
    812:         // There is no nice way to make a generic MPI type here. Really sucky
    813:         PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    814:         for(typename LabelSifter::traits::arrow_container_type::set_type::iterator ai = label._arrows.set.begin(); ai != label._arrows.set.end(); ai++, ++a) {
    815:           arrows[a*2+0] = ai->source;
    816:           arrows[a*2+1] = ai->target;
    817:         }
    818:         MPI_Send(arrows, size*2, MPIU_INT, 0, 1, label.comm());CHKERRXX(ierr);
    819:         PetscFree(arrows);CHKERRXX(ierr);
    820:       }
    821:     }
    822:     template<typename LabelSifter>
    823:     static void loadLabel(std::ifstream& fs, LabelSifter& label) {
    824:       if (label.commRank() == 0) {
    825:         // Load local
    826:         size_t numArrows;

    828:         fs >> numArrows;
    829:         for(size_t a = 0; a < numArrows; ++a) {
    830:           typename LabelSifter::traits::arrow_type::source_type source;
    831:           typename LabelSifter::traits::arrow_type::target_type target;

    833:           fs >> source;
    834:           fs >> target;
    835:           label.addArrow(typename LabelSifter::traits::arrow_type(source, target));
    836:         }
    837:         // Load and send remote
    838:         for(int p = 1; p < label.commSize(); ++p) {
    839:           PetscInt       size;
    840:           PetscInt      *arrows;

    843:           fs >> size;
    844:           MPI_Send(&size, 1, MPIU_INT, p, 1, label.comm());CHKERRXX(ierr);
    845:           PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    846:           for(PetscInt a = 0; a < size; ++a) {
    847:             fs >> arrows[a*2+0];
    848:             fs >> arrows[a*2+1];
    849:           }
    850:           MPI_Send(arrows, size*2, MPIU_INT, p, 1, label.comm());CHKERRXX(ierr);
    851:           PetscFree(arrows);CHKERRXX(ierr);
    852:         }
    853:       } else {
    854:         // Load remote
    855:         PetscInt       size;
    856:         PetscInt      *arrows;
    857:         MPI_Status     status;

    860:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, label.comm(), &status);CHKERRXX(ierr);
    861:         PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    862:         MPI_Recv(arrows, size*2, MPIU_INT, 0, 1, label.comm(), &status);CHKERRXX(ierr);
    863:         for(PetscInt a = 0; a < size; ++a) {
    864:           label.addArrow(typename LabelSifter::traits::arrow_type(arrows[a*2+0], arrows[a*2+1]));
    865:         }
    866:         PetscFree(arrows);CHKERRXX(ierr);
    867:       }
    868:     }
    869:   };
    870: } // namespace ALE

    872: #endif // ifdef included_ALE_LabelSifter_hh