Actual source code: LabelSifter.hh

  1: #ifndef included_ALE_LabelSifter_hh
  2: #define included_ALE_LabelSifter_hh

  4: #include <iostream>

  6: #ifndef  included_ALE_hh
  7: #include <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;
 27: 
 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
 64: 
 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
135: 

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;};
198: 
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:         };
216: 
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:         };
224: 
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:         };
232: 
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
254: 

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() {};
335: 
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:         };
346: 
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
353: 
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() {};
371: 
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 PetscCookie sifterType = -1;
410:       //const char        *id_name = ALE::getClassName<T>();
411:       const char        *id_name = "LabelSifter";
412:       PetscErrorCode     ierr;

414:       if (sifterType < 0) {
415:         PetscCookieRegister(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
559:     PetscErrorCode view(const char* label = NULL, bool raw = false){
561:       ostringstream txt;
563:       if(this->_debug) {
564:         std::cout << "viewing a LabelSifter, comm = " << this->comm() << ", PETSC_COMM_SELF = " << PETSC_COMM_SELF << ", commRank = " << this->commRank() << std::endl;
565:       }
566:       if(label != NULL) {
567:         PetscPrintf(this->comm(), "viewing LabelSifter: '%s'\n", label);
568:       } else {
569:         PetscPrintf(this->comm(), "viewing a LabelSifter: \n");
570:       }
571:       if(!raw) {
572:         ostringstream txt;
573:         if(this->commRank() == 0) {
574:           txt << "cap --> base:\n";
575:         }
576:         if(_arrows.set.empty()) {
577:           txt << "[" << this->commRank() << "]: empty" << std::endl;
578:         }
579:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++) {
580:           txt << "[" << this->commRank() << "]: " << ai->source << "---->" << ai->target << std::endl;
581:         }
582:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
583:         PetscSynchronizedFlush(this->comm()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
584:       }
585:       else { // if(raw)
586:         ostringstream txt;
587:         if(this->commRank() == 0) {
588:           txt << "'raw' arrow set:" << std::endl;
589:         }
590:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++)
591:         {
592:           typename traits::arrow_type arr = *ai;
593:           txt << "[" << this->commRank() << "]: " << arr << std::endl;
594:         }
595:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
596:         PetscSynchronizedFlush(this->comm()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
597:       }// if(raw)
598: 
599:       return(0);
600:     };
601:   public:
602:     //
603:     // Lattice queries
604:     //
605:     template<class targetInputSequence>
606:     Obj<typename traits::coneSequence> meet(const Obj<targetInputSequence>& targets);
607:     // unimplemented
608:     template<class sourceInputSequence>
609:     Obj<typename traits::coneSequence> join(const Obj<sourceInputSequence>& sources);
610:   public:
611:     //
612:     // Structural manipulation
613:     //
614:     void clear() {
615:       this->_arrows.set.clear();
616:     };
617:     // This is necessary to work with Completion right now
618:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q, const int dummy) {
619:       this->addArrow(p, q);
620:     };
621:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q) {
622:       this->addArrow(typename traits::arrow_type(p, q));
623:       //std::cout << "Added " << arrow_type(p, q);
624:     };
625:     virtual void addArrow(const typename traits::arrow_type& a) {
626:       this->_arrows.set.insert(a);
627:     };
628:     virtual void removeArrow(const typename traits::arrow_type& a) {
629:       // First, produce an arrow sequence for the given source, target combination.
630:       typename traits::arrowSequence::traits::index_type& arrowIndex =
631:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
632:       typename traits::arrowSequence::traits::index_type::iterator i,ii,j;
633:       i = arrowIndex.lower_bound(::boost::make_tuple(a.source,a.target));
634:       ii = arrowIndex.upper_bound(::boost::make_tuple(a.source, a.target));
635:       if (this->_debug) {
636:         std::cout << "removeArrow: attempting to remove arrow:" << a << std::endl;
637:         std::cout << "removeArrow: candidate arrows are:" << std::endl;
638:       }
639:       for(j = i; j != ii; j++) {
640:         if (this->_debug) {
641:           std::cout << " " << *j;
642:         }
643:         // Find the arrow of right color and remove it
644:         if (this->_debug) {
645:           std::cout << std::endl << "removeArrow: found:" << *j << std::endl;
646:         }
647:         arrowIndex.erase(j);
648:         break;
649:       }
650:     };

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

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

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

728:         if (!baseRestriction->getSupportSize(target) && !baseRestriction->getConeSize(target)) continue;
729:         this->addArrow(a_iter->source, target);
730:       }
731:     };
732:     int size() const {return _arrows.set.size();};
733:     int getCapSize() const {
734:       std::set<source_type> cap;
735:       for(typename traits::arrow_container_type::set_type::iterator a_iter = _arrows.set.begin(); a_iter != _arrows.set.end(); ++a_iter) {
736:         cap.insert(a_iter->source);
737:       }
738:       return cap.size();
739:     };
740:     capSequence cap() const {
741:       std::set<source_type> cap;
742:       for(typename traits::arrow_container_type::set_type::iterator a_iter = _arrows.set.begin(); a_iter != _arrows.set.end(); ++a_iter) {
743:         cap.insert(a_iter->source);
744:       }
745:       return cap;
746:     };
747:     int getBaseSize() const {
748:       std::set<target_type> base;
749:       for(typename traits::arrow_container_type::set_type::iterator a_iter = _arrows.set.begin(); a_iter != _arrows.set.end(); ++a_iter) {
750:         base.insert(a_iter->target);
751:       }
752:       return base.size();
753:     };
754:   public: // Compatibility with fixed storage variants
755:     typedef Interval<target_type> chart_type;
756:     chart_type& getChart() {static chart_type chart(0, 0); return chart;};
757:     template<typename chart_type>
758:     void setChart(const chart_type& chart) {};
759:     void setConeSize(target_type p, int s) {};
760:     void setSupportSize(source_type p, int s) {};
761:     void allocate() {};
762:     void recalculateLabel() {};
763:   }; // class LabelSifter

765:   class LabelSifterSerializer {
766:   public:
767:     template<typename LabelSifter>
768:     static void writeLabel(std::ofstream& fs, LabelSifter& label) {
769:       if (label.commRank() == 0) {
770:         // Write local
771:         fs << label._arrows.set.size() << std::endl;
772:         for(typename LabelSifter::traits::arrow_container_type::set_type::iterator ai = label._arrows.set.begin(); ai != label._arrows.set.end(); ai++) {
773:           fs << ai->source << " " << ai->target << std::endl;
774:         }
775:         // Receive and write remote
776:         for(int p = 1; p < label.commSize(); ++p) {
777:           PetscInt       size;
778:           PetscInt      *arrows;
779:           MPI_Status     status;

782:           MPI_Recv(&size, 1, MPIU_INT, p, 1, label.comm(), &status);CHKERRXX(ierr);
783:           fs << size << std::endl;
784:           PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
785:           MPI_Recv(arrows, size*2, MPIU_INT, p, 1, label.comm(), &status);CHKERRXX(ierr);
786:           for(PetscInt a = 0; a < size; ++a) {
787:             fs << arrows[a*2+0] << " " << arrows[a*2+1] << std::endl;
788:           }
789:           PetscFree(arrows);CHKERRXX(ierr);
790:         }
791:       } else {
792:         // Send remote
793:         PetscInt       size = label._arrows.set.size();
794:         PetscInt       a    = 0;
795:         PetscInt      *arrows;

798:         MPI_Send(&size, 1, MPIU_INT, 0, 1, label.comm());CHKERRXX(ierr);
799:         // There is no nice way to make a generic MPI type here. Really sucky
800:         PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
801:         for(typename LabelSifter::traits::arrow_container_type::set_type::iterator ai = label._arrows.set.begin(); ai != label._arrows.set.end(); ai++, ++a) {
802:           arrows[a*2+0] = ai->source;
803:           arrows[a*2+1] = ai->target;
804:         }
805:         MPI_Send(arrows, size*2, MPIU_INT, 0, 1, label.comm());CHKERRXX(ierr);
806:         PetscFree(arrows);CHKERRXX(ierr);
807:       }
808:     };
809:     template<typename LabelSifter>
810:     static void loadLabel(std::ifstream& fs, LabelSifter& label) {
811:       if (label.commRank() == 0) {
812:         // Load local
813:         size_t numArrows;

815:         fs >> numArrows;
816:         for(size_t a = 0; a < numArrows; ++a) {
817:           typename LabelSifter::traits::arrow_type::source_type source;
818:           typename LabelSifter::traits::arrow_type::target_type target;

820:           fs >> source;
821:           fs >> target;
822:           label.addArrow(typename LabelSifter::traits::arrow_type(source, target));
823:         }
824:         // Load and send remote
825:         for(int p = 1; p < label.commSize(); ++p) {
826:           PetscInt       size;
827:           PetscInt      *arrows;

830:           fs >> size;
831:           MPI_Send(&size, 1, MPIU_INT, p, 1, label.comm());CHKERRXX(ierr);
832:           PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
833:           for(PetscInt a = 0; a < size; ++a) {
834:             fs >> arrows[a*2+0];
835:             fs >> arrows[a*2+1];
836:           }
837:           MPI_Send(arrows, size*2, MPIU_INT, p, 1, label.comm());CHKERRXX(ierr);
838:           PetscFree(arrows);CHKERRXX(ierr);
839:         }
840:       } else {
841:         // Load remote
842:         PetscInt       size;
843:         PetscInt      *arrows;
844:         MPI_Status     status;

847:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, label.comm(), &status);CHKERRXX(ierr);
848:         PetscMalloc(size*2 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
849:         MPI_Recv(arrows, size*2, MPIU_INT, 0, 1, label.comm(), &status);CHKERRXX(ierr);
850:         for(PetscInt a = 0; a < size; ++a) {
851:           label.addArrow(typename LabelSifter::traits::arrow_type(arrows[a*2+0], arrows[a*2+1]));
852:         }
853:         PetscFree(arrows);CHKERRXX(ierr);
854:       }
855:     };
856:   };
857: } // namespace ALE

859: #endif // ifdef included_ALE_LabelSifter_hh