From 1caff921bef2789e99c24d7fc228f0af0f9fe942 Mon Sep 17 00:00:00 2001 From: Matthias Petri Date: Wed, 17 Aug 2011 12:00:07 +1000 Subject: [PATCH] - Added Intersection of k ranges over T[0..n]. (based on "New Algorithms on Wavelet Trees and Applications to Information Retrieval" by Travis Gagie, Gonzalo Navarro, Simon J. Puglisi. --- trunk/src/static/sequence/WaveletTree.cpp | 20 +++ trunk/src/static/sequence/WaveletTree.h | 2 + .../src/static/sequence/WaveletTreeNoptrs.cpp | 85 +++++++++++ trunk/src/static/sequence/WaveletTreeNoptrs.h | 50 +++++++ trunk/src/static/sequence/wt_node.h | 2 + .../src/static/sequence/wt_node_internal.cpp | 48 +++++++ trunk/src/static/sequence/wt_node_internal.h | 2 + trunk/src/static/sequence/wt_node_leaf.cpp | 8 ++ trunk/src/static/sequence/wt_node_leaf.h | 2 + trunk/tests/Makefile | 10 +- trunk/tests/testIntersection.cpp | 133 ++++++++++++++++++ 11 files changed, 359 insertions(+), 3 deletions(-) create mode 100644 trunk/tests/testIntersection.cpp diff --git a/trunk/src/static/sequence/WaveletTree.cpp b/trunk/src/static/sequence/WaveletTree.cpp index bf7ecce..0cc4621 100644 --- a/trunk/src/static/sequence/WaveletTree.cpp +++ b/trunk/src/static/sequence/WaveletTree.cpp @@ -152,4 +152,24 @@ namespace cds_static return ret; } + std::vector* WaveletTree::intersect(std::vector< pair >& ranges,size_t thres) { + std::vector* intersection = new std::vector(); + std::vector* iunmaped = new std::vector(); + + if(thres ==0) thres = ranges.size(); + + root->intersect(ranges,thres,iunmaped); + + /* unmap all symbols */ + if(iunmaped->size()>0) { + for(size_t i=0;isize();i++) { + intersection->push_back( am->unmap((*iunmaped)[i]) ); + } + } + + delete iunmaped; + return intersection; + } + }; + diff --git a/trunk/src/static/sequence/WaveletTree.h b/trunk/src/static/sequence/WaveletTree.h index 904816d..d86406e 100644 --- a/trunk/src/static/sequence/WaveletTree.h +++ b/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* intersect(std::vector< pair >& ranges,size_t threshold = 0); + protected: WaveletTree(); diff --git a/trunk/src/static/sequence/WaveletTreeNoptrs.cpp b/trunk/src/static/sequence/WaveletTreeNoptrs.cpp index 2a937c6..5d5a14a 100644 --- a/trunk/src/static/sequence/WaveletTreeNoptrs.cpp +++ b/trunk/src/static/sequence/WaveletTreeNoptrs.cpp @@ -21,6 +21,8 @@ #include +#include + 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* WaveletTreeNoptrs::intersect(std::vector< pair >& ranges,size_t thres) + { + std::vector* intersection = new std::vector(); + + if(thres == 0) { + thres = ranges.size(); + } + size_t start,end,before,left,right,newleft,newright; + uint csym,clevel; + BitSequence* bs; + + priority_queue 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 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; + } + }; diff --git a/trunk/src/static/sequence/WaveletTreeNoptrs.h b/trunk/src/static/sequence/WaveletTreeNoptrs.h index f97d54e..a449b38 100644 --- a/trunk/src/static/sequence/WaveletTreeNoptrs.h +++ b/trunk/src/static/sequence/WaveletTreeNoptrs.h @@ -34,6 +34,54 @@ using namespace std; namespace cds_static { + class intrange_t : public std::vector< pair > { + public: + intrange_t(size_t s,size_t e,uint lvl,uint sym,std::vector< pair > ranges) : + std::vector< pair >(ranges) { + this->sym = sym; + start = s; end = e; + level = lvl; + freq = 0; + for(size_t i=0;i p = ranges[i]; + freq += (p.second - p.first + 1); + } + }; + + intrange_t(size_t s,size_t e,uint lvl,uint sym) : + std::vector< pair >() { + 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(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 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* intersect(std::vector< pair >& ranges,size_t thres=0); + protected: WaveletTreeNoptrs(); diff --git a/trunk/src/static/sequence/wt_node.h b/trunk/src/static/sequence/wt_node.h index 987655c..8f29b99 100644 --- a/trunk/src/static/sequence/wt_node.h +++ b/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 >& ranges, + size_t thres,std::vector* intersection) const = 0; }; }; diff --git a/trunk/src/static/sequence/wt_node_internal.cpp b/trunk/src/static/sequence/wt_node_internal.cpp index 22a08d9..46450de 100644 --- a/trunk/src/static/sequence/wt_node_internal.cpp +++ b/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 >& ranges, + size_t thres,std::vector* 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 > ranges_zero; + std::vector< pair > ranges_one; + + for(size_t i=0;i 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(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(newleft,newright)); + } + } + left_child->intersect(ranges_zero,thres,intersection); + right_child->intersect(ranges_one,thres,intersection); + } + } }; + diff --git a/trunk/src/static/sequence/wt_node_internal.h b/trunk/src/static/sequence/wt_node_internal.h index 9d6cc92..b698a5c 100644 --- a/trunk/src/static/sequence/wt_node_internal.h +++ b/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 >& ranges, + size_t thres,std::vector* intersection) const; protected: wt_node_internal(); diff --git a/trunk/src/static/sequence/wt_node_leaf.cpp b/trunk/src/static/sequence/wt_node_leaf.cpp index d078106..390133d 100644 --- a/trunk/src/static/sequence/wt_node_leaf.cpp +++ b/trunk/src/static/sequence/wt_node_leaf.cpp @@ -81,4 +81,12 @@ namespace cds_static ret->symbol = loadValue(fp); return ret; } + + void wt_node_leaf::intersect(std::vector< pair >& ranges, + size_t thres,std::vector* intersection) const { + if(ranges.size() >= thres ) { + intersection->push_back(symbol); + } + } }; + diff --git a/trunk/src/static/sequence/wt_node_leaf.h b/trunk/src/static/sequence/wt_node_leaf.h index c5247fd..a7ebdb0 100644 --- a/trunk/src/static/sequence/wt_node_leaf.h +++ b/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 >& ranges, + size_t thres,std::vector* intersection) const; protected: wt_node_leaf(); diff --git a/trunk/tests/Makefile b/trunk/tests/Makefile index eaa7551..92db155 100644 --- a/trunk/tests/Makefile +++ b/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) diff --git a/trunk/tests/testIntersection.cpp b/trunk/tests/testIntersection.cpp new file mode 100644 index 0000000..711033b --- /dev/null +++ b/trunk/tests/testIntersection.cpp @@ -0,0 +1,133 @@ + +#include + +#include +#include +#include +#include + +using namespace std; +using namespace cds_static; + + +void testIntersectionWTNoptrs(Array & a,WaveletTreeNoptrs& wt) { + uint tests = 500; + + while(tests) { + /* create ranges */ + std::vector< pair > ranges; + uint N = 2 + rand() % 5; + for(size_t i=0;i > r2 = ranges; + + /* perform intersection */ + std::vector* isect = wt.intersect(ranges); + + /* check result */ + for(size_t j=0;jsize();j++) { + uint res = (*isect)[j]; + uint found = 0; + for(size_t i=0;i cr = r2[i]; + + for(size_t k=cr.first;k<=cr.second;k++) { + if( a[k] == res ) { + found++; + break; + } + } + } + if(found!=N) { + fprintf(stdout,"INTERSECTION ERROR!\n"); + exit(EXIT_FAILURE); + } + } + delete isect; + + tests--; + } +} + +void testIntersectionWT(Array & a,WaveletTree& wt) { + uint tests = 500; + + while(tests) { + /* create ranges */ + std::vector< pair > ranges; + uint N = 2 + rand() % 5; + for(size_t i=0;i > r2 = ranges; + + /* perform intersection */ + std::vector* isect = wt.intersect(ranges); + + /* check result */ + for(size_t j=0;jsize();j++) { + uint res = (*isect)[j]; + uint found = 0; + for(size_t i=0;i cr = r2[i]; + + for(size_t k=cr.first;k<=cr.second;k++) { + if( a[k] == res ) { + found++; + break; + } + } + } + if(found!=N) { + fprintf(stdout,"ERROR! could not find %u in enough ranges\n",res); + exit(EXIT_FAILURE); + } + } + delete isect; + + tests--; + } +} + + +int main(int argc, char ** argv) { + + srand( getpid() ); + + uint len = 250 + rand() % 50000; + uint maxv = 20 + rand() % 5; + + Array a(len,maxv); + for(uint i=0;iuse(); + mapper2->use(); + + cout << "Test 1 : Wavelet tree without pointers" << endl; + WaveletTreeNoptrs wt3(a, new BitSequenceBuilderRRR(32), mapper); + cout << "bs.size() = " << wt3.getSize() << endl; + testIntersectionWTNoptrs(a,wt3); + + cout << "Test 2 : wavelet tree with pointers" << endl; + WaveletTree wt1(a,new wt_coder_huff(a,mapper),new BitSequenceBuilderRG(20),mapper); + cout << "bs.size() = " << wt1.getSize() << endl; + testIntersectionWT(a,wt1); + mapper->unuse(); + mapper2->unuse(); + fprintf(stdout,"ALL OK\n"); + return 0; +} +