Added Intersection of k ranges over T[0..n].

trunk/src/static/sequence/WaveletTree.cpp

@@ -152,4 +152,24 @@ namespace cds_static
         return ret;
     }
 
+    std::vector<uint32_t>* WaveletTree::intersect(std::vector< pair<size_t,size_t> >& ranges,size_t thres) {
+        std::vector<uint>* intersection = new std::vector<uint>();
+        std::vector<uint>* iunmaped = new std::vector<uint>();
+
+        if(thres ==0) thres = ranges.size();
+
+        root->intersect(ranges,thres,iunmaped);
+
+        /* unmap all symbols */
+        if(iunmaped->size()>0) {
+            for(size_t i=0;i<iunmaped->size();i++) {
+                intersection->push_back( am->unmap((*iunmaped)[i]) );
+            }
+        }
+
+        delete iunmaped;
+        return intersection;
+    }
+
 };
+

trunk/src/static/sequence/WaveletTree.h

@@ -86,6 +86,8 @@ namespace cds_static
             virtual void save(ofstream & fp) const;
             static WaveletTree * load(ifstream & fp);
 
+            std::vector<uint32_t>* intersect(std::vector< pair<size_t,size_t> >& ranges,size_t threshold = 0);
+
         protected:
 
             WaveletTree();

trunk/src/static/sequence/WaveletTreeNoptrs.cpp

@@ -21,6 +21,8 @@
 
 #include <WaveletTreeNoptrs.h>
 
+#include <queue>
+
 namespace cds_static
 {
     WaveletTreeNoptrs::WaveletTreeNoptrs(const Array & a, BitSequenceBuilder * bmb, Mapper * am) : Sequence(0) {
@@ -492,4 +494,87 @@ namespace cds_static
         return occ->select1(am->map(symbol))-occ->select1(am->map(symbol)-1)+1;
     }
 
+    std::vector<uint>* WaveletTreeNoptrs::intersect(std::vector< pair<size_t,size_t> >& ranges,size_t thres)
+    {
+        std::vector<uint>* intersection = new std::vector<uint>();
+
+        if(thres == 0) {
+            thres = ranges.size();
+        }
+        size_t start,end,before,left,right,newleft,newright;
+        uint csym,clevel;
+        BitSequence* bs;
+
+        priority_queue<intrange_t> heap;
+        heap.push( intrange_t(0,n-1,0,0, ranges) );
+
+        while( ! heap.empty() ) {
+            intrange_t cur_ranges = heap.top();
+            heap.pop();
+
+            /* check if we have non-empty valid ranges */
+            if( cur_ranges.contains_empty() ) continue;
+
+            csym = cur_ranges.sym;
+            clevel = cur_ranges.level;
+            if( clevel == height ) {
+                 if( thres <= cur_ranges.size() ) {
+                     /* we found a symbol */
+                     intersection->push_back( am->unmap(csym) );
+                 }
+            } else {
+                /* map each range to the corresponding range at level + 1 */
+                start = cur_ranges.start;
+                end = cur_ranges.end;
+                bs = bitstring[clevel];
+
+                if(start == 0) before = 0;
+                else before = bs->rank1(start-1);
+
+                /* new interval bounds */
+                size_t newintbound = end - (bs->rank1(end)-before);
+
+                intrange_t ranges_zero(start,newintbound,clevel+1,csym);
+                intrange_t ranges_one(newintbound+1,end,clevel+1, set(csym,clevel) );
+
+                for(size_t i=0;i<cur_ranges.size();i++) {
+                    pair<size_t,size_t> r = cur_ranges[i];
+                    left = r.first;
+                    right = r.second;
+
+                    /* number of 1s before T[l..r] */
+                    size_t rank_before_left = bs->rank1(start+left-1);
+                    /* number of 1s before T[r] */
+                    size_t rank_before_right = bs->rank1(start+right);
+                    /* number of 1s in T[l..r] */
+                    size_t num_ones = rank_before_right - rank_before_left;
+                    /* number of 1s in T[l..r] */
+                    size_t num_zeros = (right-left+1) - num_ones;
+
+                    if(num_ones) {
+                        /* number of 1s before T[l..r] within the current node */
+                        newleft = rank_before_left - before;
+                        /* number of 1s in T[l..r] */
+                        newright = rank_before_right - before - 1;
+
+                        ranges_one.addRange(newleft,newright);
+                    }
+
+                    if(num_zeros) {
+                        /* number of zeros before T[l..r] within the current 
+                         * node */
+                        newleft = left - (rank_before_left - before);
+                        /* number of zeros in T[l..r] + left bound */
+                        newright = right - (rank_before_right - before);
+
+                        ranges_zero.addRange(newleft,newright);
+                    }
+                }
+                heap.push(ranges_zero);
+                heap.push(ranges_one);
+            }
+        }
+        return intersection;
+    }
+
 };

trunk/src/static/sequence/WaveletTreeNoptrs.h

@@ -34,6 +34,54 @@ using namespace std;
 
 namespace cds_static
 {
+    class intrange_t : public std::vector< pair<size_t,size_t> > {
+        public:
+            intrange_t(size_t s,size_t e,uint lvl,uint sym,std::vector< pair<size_t,size_t> > ranges) :
+                    std::vector< pair<size_t,size_t> >(ranges) {
+                this->sym = sym;
+                start = s; end = e;
+                level = lvl;
+                freq = 0;
+                for(size_t i=0;i<ranges.size();i++) {
+                    pair<size_t,size_t> p = ranges[i];
+                    freq += (p.second - p.first + 1);
+                }
+            };
+
+            intrange_t(size_t s,size_t e,uint lvl,uint sym) :
+                    std::vector< pair<size_t,size_t> >() {
+                this->sym = sym;
+                start = s; end = e;
+                level = lvl;
+                freq = 0;
+            };
+
+            void addRange( size_t sp, size_t ep ) {
+                freq += (ep-sp+1);
+                push_back( pair<size_t,size_t>(sp,ep) );
+            };
+
+            bool operator<(const intrange_t& r) const
+            {
+                return ( (freq) < (r.freq));
+            };
+
+            bool contains_empty() {
+                if( size() == 0 ) return true;
+                for(size_t i=0;i<size();i++) {
+                    pair<size_t,size_t> interval = (*this)[i];
+                    if( interval.first > interval.second ) return true;
+                }
+                return false;
+            };
+        public:
+            uint sym;
+            size_t start;
+            size_t end;
+            uint level;
+            size_t freq;
+    };
+
 
     class WaveletTreeNoptrs : public Sequence
     {
@@ -74,6 +122,8 @@ namespace cds_static
             virtual void save(ofstream & fp) const;
             static WaveletTreeNoptrs * load(ifstream & fp);
 
+            std::vector<uint32_t>* intersect(std::vector< pair<size_t,size_t> >& ranges,size_t thres=0);
+
         protected:
             WaveletTreeNoptrs();
 

trunk/src/static/sequence/wt_node.h

@@ -50,6 +50,8 @@ namespace cds_static
             virtual size_t getSize() const = 0;
             virtual void save(ofstream & fp) const = 0;
             static wt_node * load(ifstream & fp);
+            virtual void intersect(std::vector< pair<size_t,size_t> >& ranges,
+                    size_t thres,std::vector<uint>* intersection) const = 0;
     };
 
 };

trunk/src/static/sequence/wt_node_internal.cpp

@@ -260,4 +260,52 @@ namespace cds_static
         return ret;
     }
 
+    void wt_node_internal::intersect(std::vector< pair<size_t,size_t> >& ranges,
+            size_t thres,std::vector<uint>* intersection) const {
+        if(ranges.size() < thres) {
+            return;
+        } else {
+            /* map each range to the corresponding range at level + 1 */
+            size_t left,right,newleft,newright;
+
+            std::vector< pair<size_t,size_t> > ranges_zero;
+            std::vector< pair<size_t,size_t> > ranges_one;
+
+            for(size_t i=0;i<ranges.size();i++) {
+                pair<size_t,size_t> r = ranges[i];
+                left = r.first;
+                right = r.second;
+
+                /* number of 1s before T[l..r] */
+                size_t rank_before_left = bitmap->rank1(left-1);
+                /* number of 1s before T[r] */
+                size_t rank_before_right = bitmap->rank1(right);
+                /* number of 1s in T[l..r] */
+                size_t num_ones = rank_before_right - rank_before_left;
+                /* number of 1s in T[l..r] */
+                size_t num_zeros = (right-left+1) - num_ones;
+
+                if(num_ones) {
+                    /* number of 1s before T[l..r] within the current node */
+                    newleft = rank_before_left;
+                    /* number of 1s in T[l..r] */
+                    newright = newleft+num_ones-1;
+
+                    ranges_one.push_back( pair<size_t,size_t>(newleft,newright) );
+                }
+
+                if(num_zeros) {
+                    /* number of zeros before T[l..r] within the current node */
+                    newleft = left - rank_before_left;
+                    /* number of zeros in T[l..r] */
+                    newright = newleft+num_zeros-1;
+
+                    ranges_zero.push_back( pair<size_t,size_t>(newleft,newright));
+                }
+            }
+            left_child->intersect(ranges_zero,thres,intersection);
+            right_child->intersect(ranges_one,thres,intersection);
+        }
+    }
 };
+

trunk/src/static/sequence/wt_node_internal.h

@@ -49,6 +49,8 @@ namespace cds_static
             virtual size_t getSize() const;
             virtual void save(ofstream & fp) const;
             static wt_node_internal * load(ifstream & fp);
+            virtual void intersect(std::vector< pair<size_t,size_t> >& ranges,
+                    size_t thres,std::vector<uint>* intersection) const;
 
         protected:
             wt_node_internal();

trunk/src/static/sequence/wt_node_leaf.cpp

@@ -81,4 +81,12 @@ namespace cds_static
         ret->symbol = loadValue<uint>(fp);
         return ret;
     }
+
+    void wt_node_leaf::intersect(std::vector< pair<size_t,size_t> >& ranges,
+            size_t thres,std::vector<uint>* intersection) const {
+        if(ranges.size() >= thres ) {
+            intersection->push_back(symbol);
+        }
+    }
 };
+

trunk/src/static/sequence/wt_node_leaf.h

@@ -47,6 +47,8 @@ namespace cds_static
             virtual size_t getSize() const;
             virtual void save(ofstream & fp) const;
             static wt_node_leaf * load(ifstream & fp);
+            virtual void intersect(std::vector< pair<size_t,size_t> >& ranges,
+                    size_t thres,std::vector<uint>* intersection) const;
 
         protected:
             wt_node_leaf();

trunk/tests/Makefile

@@ -2,12 +2,12 @@
 CPP=g++
 
 #CPPFLAGS=-g3 -Wall -O0 -I../includes/ 
-#CPPFLAGS=-O3 -Wall -DNDEBUG -I../includes/
+CPPFLAGS=-O3 -Wall -DNDEBUG -I../includes/
 INCS=-I../includes/
 LIB=../lib/libcds.a
 
-OBJECTS= testArray.o testBitSequence.o testSequence.o testHuffman.o testTextIndex.o testLCP.o testNPR.o testSuffixTree.o timeSequence.o toArray.o toArray2.o
-BIN= testArray testBitSequence testSequence testHuffman testTextIndex testLCP testNPR testSuffixTree timeSequence toArray toArray2
+OBJECTS= testArray.o testBitSequence.o testSequence.o testHuffman.o testTextIndex.o testLCP.o testNPR.o testSuffixTree.o timeSequence.o toArray.o toArray2.o testIntersection.o
+BIN= testArray testBitSequence testSequence testHuffman testTextIndex testLCP testNPR testSuffixTree timeSequence toArray toArray2 testIntersection
 
 %.o: %.cpp
 	@echo " [C++] Compiling $<"
@@ -29,6 +29,10 @@ testBitSequence:
 	@echo " [LNK] Building testBitSequence"
 	@$(CPP) $(CPPFLAGS) -o testBitSequence testBitSequence.o $(LIB)
 
+testIntersection:
+	@echo " [LNK] Building testIntersection"
+	@$(CPP) $(CPPFLAGS) -o testIntersection testIntersection.o $(LIB)
+
 testArray:
 	@echo " [LNK] Building testArray"
 	@$(CPP) $(CPPFLAGS) -o testArray testArray.o $(LIB)