{
makewavelet(bwt); // Deletes bwt!
bwt = 0;
+
+ CSA::DeltaVector::Iterator notIndexedIt(notIndexed);
// Make sampling tables
- maketables(numberOfSamples_, tsType, notIndexed, niText);
+ maketables(numberOfSamples_, tsType, notIndexedIt, niText);
}
bool FMIndex::EmptyText(DocId k) const
#endif
}
-void FMIndex::maketables(ulong sampleLength, char tsType, CSA::DeltaVector & notIndexed, const string & niText)
+void FMIndex::maketables(ulong sampleLength, char tsType, CSA::DeltaVector::Iterator & notIndexedIt, const string & niText)
{
// Calculate BWT end-marker position (of last inserted text)
{
if (prevTextId == numberOfTexts)
{
prevTextId = 0;
- while (notIndexed.isSet(prevTextId))
+ while (notIndexedIt.isSet(prevTextId))
++ prevTextId;
// Now prevTextId points to the first indexed Doc ID.
}
/**
* Insert non-indexed texts
*/
- while (notIndexed.isSet(textId))
+ while (notIndexedIt.isSet(textId))
{
do {
tsbuilder[tsb_i] = *nit_i;
// LF-mapping from '\0' does not work with this (pseudo) BWT.
// Correct LF-mapping to the last char of the previous text:
- p = textId - notIndexed.rank(textId);
+ p = textId - notIndexedIt.rank(textId);
}
else // Now c != '\0', do LF-mapping:
p = C[c]+alphabetrank_i_tmp-1;
}
- while (textId > 0 && notIndexed.isSet(textId-1))
+ while (textId > 0 && notIndexedIt.isSet(textId-1))
{
do {
-- tsb_i;
for(ulong i = 0; i < sampleLength; i ++) {
// Find next sampled text position
while ((posOfSuccEndmarker - x) % samplerate != 0
- || notIndexed.isSet(textId)) // Loop over non-indexed
+ || notIndexedIt.isSet(textId)) // Loop over non-indexed
{
--x;
assert(x != ~0lu);
// Following methods are not part of the public API
uchar * BWT(uchar *);
void makewavelet(uchar *);
- void maketables(ulong, char, CSA::DeltaVector &, const std::string &);
+ void maketables(ulong, char, CSA::DeltaVector::Iterator &, const std::string &);
DocId DocIdAtTextPos(BlockArray*, TextPosition) const;
ulong Search(uchar const *, TextPosition, TextPosition *, TextPosition *) const;
ulong Search(uchar const *, TextPosition, TextPosition *, TextPosition *, DocId, DocId) const;
}
p_->notIndexed->setBit(p_->numberOfTexts); // FIXME CSA::DeltaVector can not be all 0's
- CSA::DeltaVector deltav = CSA::DeltaVector(*p_->notIndexed, p_->numberOfTexts+1);
+ CSA::DeltaVector deltav(*p_->notIndexed, p_->numberOfTexts+1);
delete p_->notIndexed;
p_->notIndexed = 0;
--- /dev/null
+#include "incbwt/rlcsa_builder.h"
+#include "RLCSABuilder.h"
+#include "RLCSAWrapper.h"
+
+namespace SXSI
+{
+
+// Using pimpl idiom to hide RLCSABuilder*
+struct TCBuilderRep
+{
+ unsigned samplerate;
+ CSA::RLCSABuilder * sa;
+
+ ulong n;
+ unsigned numberOfTexts;
+ bool insertAllowed;
+};
+
+/**
+ * Init text collection
+ */
+RLCSABuilder::RLCSABuilder(unsigned samplerate, ulong estimatedInputLength)
+ : p_(new struct TCBuilderRep())
+{
+ p_->n = 0;
+ p_->samplerate = samplerate;
+ if (samplerate == 0)
+ p_->samplerate = TEXTCOLLECTION_DEFAULT_SAMPLERATE;
+
+ p_->numberOfTexts = 0;
+ p_->insertAllowed = true;
+
+ CSA::usint rlcsa_block_size = CSA::RLCSA_BLOCK_SIZE.second;
+ CSA::usint rlcsa_sample_rate = p_->samplerate;
+
+ // Parameters for RLCSA: 32 bytes, samples, buffer size n/10 bytes.
+ // Buffer size is always at least 15MB:
+ if (estimatedInputLength < TEXTCOLLECTION_DEFAULT_INPUT_LENGTH)
+ estimatedInputLength = TEXTCOLLECTION_DEFAULT_INPUT_LENGTH;
+ p_->sa = new CSA::RLCSABuilder(rlcsa_block_size, rlcsa_sample_rate, estimatedInputLength/10);
+ assert(p_->sa->isOk());
+}
+
+RLCSABuilder::~RLCSABuilder()
+{
+ delete p_->sa;
+ delete p_;
+}
+
+void RLCSABuilder::InsertText(uchar const * text, bool index)
+{
+ assert(index);
+ if (!index)
+ {
+ std::cerr << "SWCSABuilder::InsertText(): The implementation of SWCSA does not support non-indexed texts"
+ << std::endl << "Use the default (FMIndex) text collection instead." << std::endl;
+ std::exit(1);
+ }
+
+ if (!p_->insertAllowed)
+ {
+ std::cerr << "RLCSABuilder::InsertText() error: new text can not be inserted after InitTextCollection() call!" << std::endl;
+ std::exit(1);
+ }
+
+ TextCollection::TextPosition m = std::strlen((char *)text) + 1;
+ if (m > 1)
+ {
+ p_->n += m;
+ p_->numberOfTexts ++;
+
+ p_->sa->insertSequence((char*)text, m-1, 0);
+ assert(p_->sa->isOk());
+ }
+ else
+ {
+ // FIXME indexing empty texts
+ std::cerr << "RLCSABuilder::InsertText() error: can not index empty texts!" << std::endl;
+ exit(1);
+ }
+}
+
+TextCollection * RLCSABuilder::InitTextCollection(char type)
+{
+ p_->insertAllowed = false; // Disable future insertions
+ assert(type == TextStorage::TYPE_PLAIN_TEXT);
+ if (type != TextStorage::TYPE_PLAIN_TEXT)
+ {
+ std::cerr << "RLCSABuilder::InitTextCollection(): The implementation of RLCSA supports only TextStorage::TYPE_PLAIN_TEXT"
+ << std::endl << "Use the default (FMIndex) text collection instead." << std::endl;
+ std::exit(1);
+ }
+
+ TextCollection *result = new RLCSAWrapper(p_->sa->getRLCSA());
+ delete p_->sa;
+ p_->sa = 0;
+ return result;
+}
+}
--- /dev/null
+/******************************************************************************
+ * Copyright (C) 2009 by Niko Valimaki <nvalimak@cs.helsinki.fi> *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU Lesser General Public License as published *
+ * by the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU Lesser General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU Lesser General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ******************************************************************************/
+
+#ifndef _SXSI_RLCSABuilder_h_
+#define _SXSI_RLCSABuilder_h_
+
+#include "TextCollectionBuilder.h"
+#include "TextStorage.h"
+#include "Tools.h" // Defines ulong and uchar.
+
+#include <string>
+#include <vector>
+#include <utility> // Defines std::pair.
+#include <cstring> // Defines std::strlen, added by Kim
+
+namespace SXSI
+{
+ struct TCBuilderRep; // Pimpl
+
+ /**
+ * Build an instance of the TextCollection class.
+ */
+ class RLCSABuilder : public TextCollectionBuilder
+ {
+ public:
+ RLCSABuilder(unsigned samplerate, ulong estimatedInputLength);
+
+ virtual ~RLCSABuilder();
+
+ /**
+ * Insert text
+ *
+ * Must be a zero-terminated string from alphabet [1,255].
+ * Can not be called after makeStatic().
+ * The i'th text insertion gets an identifier value i-1.
+ * In other words, document identifiers start from 0.
+ *
+ * Second parameter tells if the text will be added to the
+ * index also. If false, text is added only to the TextCollection
+ * and can not be searched for.
+ */
+ virtual void InsertText(uchar const *, bool index = true);
+
+ /**
+ * Make static
+ *
+ * Convert to a static collection.
+ * New texts can not be inserted after this operation.
+ *
+ * TextStorage type defaults to TYPE_PLAIN_TEXT, another
+ * possible type is TYPE_LZ_INDEX.
+ */
+ virtual TextCollection * InitTextCollection(char type = TextStorage::TYPE_PLAIN_TEXT);
+
+ private:
+ RLCSABuilder();
+
+ // Using Pimpl idiom to hide RLCSA implementation.
+ struct TCBuilderRep * p_;
+
+ // No copy constructor or assignment
+ RLCSABuilder(RLCSABuilder const&);
+ RLCSABuilder& operator = (RLCSABuilder const&);
+ };
+}
+#endif
--- /dev/null
+/******************************************************************************
+ * Copyright (C) 2010 by Niko Välimäki *
+ * *
+ * RLCSA implementation for the TextCollection interface *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU Lesser General Public License as published *
+ * by the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU Lesser General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU Lesser General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. *
+ *****************************************************************************/
+
+#ifndef _RLCSAWrapper_H_
+#define _RLCSAWrapper_H_
+
+#include "TextCollection.h"
+
+#include "incbwt/rlcsa.h"
+
+// Re-define word size to ulong:
+#undef W
+#if __WORDSIZE == 64
+# define W 64
+#else
+# define W 32
+#endif
+
+#include <set>
+#include <string>
+
+namespace SXSI
+{
+
+/**
+ * Partial implementation of the TextCollection interface
+ *
+ * Supports index construction, save, load and simple search.
+ * Use FMIndex implementation for full support.
+ */
+class RLCSAWrapper : public SXSI::TextCollection
+{
+public:
+ RLCSAWrapper(const CSA::RLCSA* index)
+ : rlcsa(index)
+ { /* NOP */ }
+
+ ~RLCSAWrapper()
+ {
+ delete rlcsa; rlcsa = 0;
+ }
+
+ bool EmptyText(DocId k) const
+ {
+ return false; // Empty texts are not indexed
+ }
+
+ /**
+ * Extracting one text.
+ *
+ * Call DeleteText() for each pointer returned by GetText()
+ * to avoid possible memory leaks.
+ */
+ uchar * GetText(DocId i) const
+ {
+ return rlcsa->display(i);
+ }
+ void DeleteText(uchar *text) const
+ {
+ delete [] text;
+ }
+
+ /**
+ * Returns a pointer to the beginning of texts i, i+1, ..., j.
+ * Texts are separated by a '\0' byte.
+ *
+ * Call DeleteText() for each pointer returned by GetText()
+ * to avoid possible memory leaks.
+ */
+ uchar * GetText(DocId i, DocId j) const
+ {
+ std::cerr << "RLCSAWrapper::GetText(i,j): unsupported method!" << std::endl
+ << "Use the default (FMIndex) text collection instead!" << std::endl;
+ std::exit(1);
+ }
+
+ bool IsPrefix(uchar const *) const { unsupported(); return false; };
+ bool IsSuffix(uchar const *) const { unsupported(); return false; };
+ bool IsEqual(uchar const *) const { unsupported(); return false; };
+ bool IsContains(uchar const *) const { unsupported(); return false; };
+ bool IsLessThan(uchar const *) const { unsupported(); return false; };
+
+ bool IsPrefix(uchar const *, DocId, DocId) const { unsupported(); return false; };
+ bool IsSuffix(uchar const *, DocId, DocId) const { unsupported(); return false; };
+ bool IsEqual(uchar const *, DocId, DocId) const { unsupported(); return false; };
+ bool IsContains(uchar const *, DocId, DocId) const { unsupported(); return false; };
+ bool IsLessThan(uchar const *, DocId, DocId) const { unsupported(); return false; };
+
+ ulong Count(uchar const *pattern) const
+ {
+// std::pair<TextPosition, TextPosition> p = backwards(pattern);
+ CSA::pair_type p = rlcsa->count(std::string((char const *)pattern));
+
+ if (p.first <= p.second)
+ return p.second - p.first + 1;
+ else
+ return 0;
+ }
+
+ unsigned CountPrefix(uchar const *) const { unsupported(); return 0; };
+ unsigned CountSuffix(uchar const *) const { unsupported(); return 0; };
+ unsigned CountEqual(uchar const *) const { unsupported(); return 0; };
+ unsigned CountContains(uchar const *) const { unsupported(); return 0; };
+ unsigned CountLessThan(const unsigned char*) const { unsupported(); return 0; };
+
+ unsigned CountPrefix(uchar const *, DocId, DocId) const { unsupported(); return 0; };
+ unsigned CountSuffix(uchar const *, DocId, DocId) const { unsupported(); return 0; };
+ unsigned CountEqual(uchar const *, DocId, DocId) const { unsupported(); return 0; };
+ unsigned CountContains(uchar const *, DocId, DocId) const { unsupported(); return 0; };
+ unsigned CountLessThan(uchar const *, DocId, DocId) const { unsupported(); return 0; };
+
+ // Definition of document_result is inherited from SXSI::TextCollection.
+ document_result Prefix(uchar const *) const { unsupported(); return document_result(); };
+ document_result Suffix(uchar const *) const { unsupported(); return document_result(); };
+ document_result Equal(uchar const *) const { unsupported(); return document_result(); };
+ document_result Contains(uchar const *pattern) const
+ {
+ //std::pair<TextPosition,TextPosition> p = backwards(pattern);
+ CSA::pair_type p = rlcsa->count(std::string((char const *)pattern));
+ if (p.first > p.second)
+ return document_result();
+
+ document_result dr;
+ dr.reserve(p.second - p.first + 2);
+ for (ulong i = p.first; i <= p.second; ++i)
+ dr.push_back(rlcsa->getSequenceForPosition(rlcsa->locate(i)));
+ return dr;
+ }
+ document_result LessThan(uchar const *) const { unsupported(); return document_result(); };
+ document_result KMismaches(uchar const *, unsigned) const { unsupported(); return document_result(); };
+ document_result KErrors(uchar const *, unsigned) const { unsupported(); return document_result(); };
+
+ document_result Prefix(uchar const *, DocId, DocId) const { unsupported(); return document_result(); };
+ document_result Suffix(uchar const *, DocId, DocId) const { unsupported(); return document_result(); };
+ document_result Equal(uchar const *, DocId, DocId) const { unsupported(); return document_result(); };
+ document_result Contains(uchar const *, DocId, DocId) const { unsupported(); return document_result(); };
+ document_result LessThan(uchar const *, DocId, DocId) const { unsupported(); return document_result(); };
+
+ // Definition of full_result is inherited from SXSI::TextCollection.
+ full_result FullContains(uchar const *) const { unsupported(); return full_result(); };
+ full_result FullContains(uchar const *, DocId, DocId) const { unsupported(); return full_result(); };
+ full_result FullKMismatches(uchar const *, unsigned) const { unsupported(); return full_result(); };
+ full_result FullKErrors(uchar const *, unsigned) const { unsupported(); return full_result(); };
+
+ // Index from/to disk
+ RLCSAWrapper(FILE *file, char const *filename)
+ : rlcsa(new CSA::RLCSA(std::string(filename)))
+ { /* NOP */ }
+
+ void Save(FILE *file, char const *filename) const
+ {
+ const char type = 'R';
+ // Saving type info:
+ if (std::fwrite(&type, 1, 1, file) != 1)
+ throw std::runtime_error("RLCSAWrapper::Save(): file write error (type flag).");
+
+ this->rlcsa->writeTo(std::string(filename));
+ }
+
+private:
+ const CSA::RLCSA* rlcsa;
+
+/* std::pair<TextPosition, TextPosition> backwards(uchar const *pattern) const
+ {
+ TextPosition i = std::strlen((char const *)pattern) - 1;
+ int c = (int)pattern[i];
+
+ std::cerr << "\nFirst symb = " << (char)c << std::endl;
+
+ TextPosition sp = rlcsa->C(c);
+ TextPosition ep = rlcsa->C(c+1)-1;
+ printf("i = %lu, c = %c, sp = %lu, ep = %lu\n", i, pattern[i], sp, ep);
+ while (sp<=ep && i>=1)
+ {
+ c = (int)pattern[--i];
+ sp = LF(c, sp);
+ ep = LF(c, ep+1)-1;
+ printf("new: c = %c, sp = %lu, ep = %lu\n", pattern[i], sp, ep);
+ }
+
+ uchar* found = rlcsa->display(sp, std::strlen((char const *)pattern), 5, i);
+ std::cerr << "found: " << found << " (" << i << std::endl;
+
+ return std::make_pair(sp, ep);
+ }*/
+
+ inline TextCollection::TextPosition
+ LF(uchar c, TextPosition i) const
+ {
+ if(i == (TextPosition)-1 || i < this->rlcsa->getNumberOfSequences())
+ { return this->rlcsa->C(c); }
+ return this->rlcsa->LF(i - this->rlcsa->getNumberOfSequences(), c);
+ }
+
+
+ void unsupported() const
+ {
+ std::cerr << std::endl << "-------------------------------------------------------------\n"
+ << "RLCSAWrapper: unsupported method!\nSee RLCSAWrapper.h for more details.\n"
+ << "The default index (FMIndex) implements this method!" << std::endl;
+ std::exit(5);
+ }
+}; // class RLCSAWrapper
+
+} // namespace SXSI
+
+#endif
/******************************************************************************
* Copyright (C) 2010 by Niko Välimäki *
* *
- * FMIndex implementation for the TextCollection interface *
+ * SWCSA implementation for the TextCollection interface *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
public:
SWCSAWrapper(uchar * text, ulong length, unsigned samplerate,
unsigned numberOfTexts_)
- : index(0), offsets(0), n(length), seSize(0), numberOfTexts(numberOfTexts_)
+ : index(0), offsets(0), offit(0), n(length), seSize(0), numberOfTexts(numberOfTexts_)
{
// Inicializes the arrays used to detect if a char is valid or not.
StartValid();
fprintf(stderr,"\n-----------------------\nnumber of words = %lu\n------------------------\n", seSize);
encoder.setBit(seSize);
offsets = new CSA::DeltaVector(encoder, seSize+1);
- std::cerr << "Number of texts: " << numberOfTexts << " vs " << offsets->rank(seSize - 1) << std::endl;
+ offit = new CSA::DeltaVector::Iterator(*offsets);
+
+ std::cerr << "Number of texts: " << numberOfTexts << " vs " << offit->rank(seSize - 1) << std::endl;
char opt[100];
snprintf(opt, 99, "sA=%u;sAinv=%u;sPsi=%u", samplerate, samplerate, samplerate);
std::exit(r);
}
index = 0;
+ delete offit; offit = 0;
delete offsets; offsets = 0;
}
{
ulong from, to, l;
uchar *text;
- from = offsets->select(i);
- to = offsets->select(i+1); // ADD one 1-bit in to end!!!
+ from = offit->select(i);
+ to = offit->select(i+1); // ADD one 1-bit in to end!!!
int r = extractWords(index, from, to, &text, &l);
if (r)
document_result dr;
dr.reserve(numocc+1);
for (ulong i = 0; i < numocc; ++i)
- dr.push_back(offsets->rank(occ[i])-1);
+ dr.push_back(offit->rank(occ[i])-1);
free(occ);
return dr;
// Index from/to disk
SWCSAWrapper(FILE *file, char const *filename)
- : index(0), offsets(0), n(0), seSize(0), numberOfTexts(0)
+ : index(0), offsets(0), offit(0), n(0), seSize(0), numberOfTexts(0)
{
uchar verFlag = 0;
if (std::fread(&verFlag, 1, 1, file) != 1)
if (std::fread(&(this->numberOfTexts), sizeof(unsigned), 1, file) != 1)
throw std::runtime_error("FMIndex::Load(): file read error (numberOfTexts).");
- offsets = new CSA::DeltaVector(file);
+ offsets = new CSA::DeltaVector(file);
+ offit = new CSA::DeltaVector::Iterator(*offsets);
// FIXME Const correctness is broken!
int r = load_index((char *)filename, &index);
private:
void *index;
CSA::DeltaVector *offsets;
+ CSA::DeltaVector::Iterator *offit;
TextPosition n;
ulong seSize;
#include "TextCollection.h"
#include "FMIndex.h"
#include "SWCSAWrapper.h"
+#include "RLCSAWrapper.h"
namespace SXSI
{
case 'W':
return new SWCSAWrapper(fp, filename);
break;
+ case 'R':
+ return new RLCSAWrapper(fp, filename);
+ break;
}
std::cerr << "TextCollection::Load(): invalid save file version or corrupted input file." << std::endl;
#include "TextCollectionBuilder.h"
#include "FMIndexBuilder.h"
#include "SWCSABuilder.h"
+#include "RLCSABuilder.h"
namespace SXSI
{
case index_type_swcsa:
return new SWCSABuilder(samplerate);
break;
+ case index_type_rlcsa:
+ return new RLCSABuilder(samplerate, estimatedInputLength);
+ break;
}
std::cerr << "TextCollectionBuilder::create(): unknown type given: expecting enum value, type = " << type << std::endl;
std::exit(2);
// Index type defaults to FM-index.
// SWCSA can be used for natural language inputs.
// NB: Current SWCSA uses a lot of memory during construction!
- enum index_type_t { index_type_default, index_type_swcsa };
+ enum index_type_t { index_type_default, index_type_swcsa, index_type_rlcsa };
static TextCollectionBuilder* create(unsigned samplerate = TEXTCOLLECTION_DEFAULT_SAMPLERATE,
index_type_t type = index_type_default,
}
TextStorage::TextStorage(std::FILE * file)
- : n_(0), offsets_(0), numberOfTexts_(0)
+ : n_(0), offsets_(0), offit_(0), numberOfTexts_(0)
{
if (std::fread(&(this->n_), sizeof(TextPosition), 1, file) != 1)
throw std::runtime_error("TextStorage::Load(): file read error (n_).");
throw std::runtime_error("TextStorage::Load(): file read error (numberOfTexts_).");
offsets_ = new CSA::DeltaVector(file);
+ offit_ = new CSA::DeltaVector::Iterator(*offsets_);
}
void TextStorage::Save(FILE *file, char type) const
{
assert(docId < (TextCollection::DocId)numberOfTexts_);
- TextPosition from = offsets_->select(docId);
+ TextPosition from = offit_->select(docId);
TextPosition to = 0;
if (docId < (TextCollection::DocId)numberOfTexts_ - 1)
- to = offsets_->select(docId + 1) - 1;
+ to = offit_->select(docId + 1) - 1;
else
to = n_-1;
assert(i < (TextCollection::DocId)numberOfTexts_);
assert(j < (TextCollection::DocId)numberOfTexts_);
- TextPosition from = offsets_->select(i);
+ TextPosition from = offit_->select(i);
TextPosition to = 0;
if (j < (TextCollection::DocId)numberOfTexts_ - 1)
- to = offsets_->select(j + 1) - 1;
+ to = offit_->select(j + 1) - 1;
else
to = n_-1;
while (i < j && i < (TextCollection::DocId)numberOfTexts_)
{
++i;
- text[offsets_->select(i) - 1 - from] = 0;
+ text[offit_->select(i) - 1 - from] = 0;
}
text[l-1] = 0;
return text;
virtual ~TextStorage()
{
+ delete offit_;
+ offit_ = 0;
delete offsets_;
- offsets_ = 0;
+ offsets_ = 0;
}
TextCollection::DocId DocIdAtTextPos(TextCollection::TextPosition i) const
{
assert(i < n_);
- return offsets_->rank(i)-1;
+ return offit_->rank(i)-1;
}
TextCollection::TextPosition TextStartPos(TextCollection::DocId i) const
{
assert(i < (TextCollection::DocId)numberOfTexts_);
- return offsets_->select(i);
+ return offit_->select(i);
}
bool IsEndmarker(TextCollection::TextPosition i) const
assert(i < n_);
if (i >= n_ - 1)
return true;
- return offsets_->isSet(i+1);
+ return offit_->isSet(i+1);
}
const static CSA::usint DV_BLOCK_SIZE = 32;
TextStorage(uchar const * text, TextPosition n)
- : n_(n), offsets_(0), numberOfTexts_(0)
+ : n_(n), offsets_(0), offit_(0), numberOfTexts_(0)
{
// Delta encoded bitvector of text offsets.
CSA::DeltaEncoder encoder(DV_BLOCK_SIZE);
offsets_ = new CSA::DeltaVector(encoder, n_);
- numberOfTexts_ = offsets_->rank(n_ - 1);
+ offit_ = new CSA::DeltaVector::Iterator(*(offsets_));
+ numberOfTexts_ = offit_->rank(n_ - 1);
}
TextStorage(std::FILE *);
TextPosition n_;
CSA::DeltaVector *offsets_;
+ CSA::DeltaVector::Iterator *offit_;
TextPosition numberOfTexts_;
};
{
assert(docId < (TextCollection::DocId)numberOfTexts_);
- TextPosition offset = offsets_->select(docId);
+ TextPosition offset = offit_->select(docId);
return &text_[offset];
}
assert(i < (TextCollection::DocId)numberOfTexts_);
assert(j < (TextCollection::DocId)numberOfTexts_);
- TextPosition offset = offsets_->select(i);
+ TextPosition offset = offit_->select(i);
return &text_[offset];
}
incbwt/bits/bitvector.h incbwt/bits/../misc/definitions.h \
incbwt/bits/bitbuffer.h incbwt/bits/rlevector.h incbwt/sasamples.h \
incbwt/misc/definitions.h incbwt/bits/bitbuffer.h \
- incbwt/bits/deltavector.h incbwt/misc/parameters.h \
- incbwt/misc/definitions.h incbwt/bits/deltavector.h FMIndexBuilder.h \
- TextCollectionBuilder.h TextCollection.h Tools.h TextStorage.h FMIndex.h \
- BitRank.h BlockArray.h ArrayDoc.h
+ incbwt/bits/deltavector.h incbwt/lcpsamples.h incbwt/bits/array.h \
+ incbwt/misc/utils.h incbwt/misc/definitions.h incbwt/misc/parameters.h \
+ incbwt/bits/deltavector.h FMIndexBuilder.h TextCollectionBuilder.h \
+ TextCollection.h Tools.h TextStorage.h FMIndex.h BitRank.h BlockArray.h \
+ ArrayDoc.h
HeapProfiler.o: HeapProfiler.cpp HeapProfiler.h
+RLCSABuilder.o: RLCSABuilder.cpp incbwt/rlcsa_builder.h incbwt/rlcsa.h \
+ incbwt/bits/vectors.h incbwt/bits/deltavector.h incbwt/bits/bitvector.h \
+ incbwt/bits/../misc/definitions.h incbwt/bits/bitbuffer.h \
+ incbwt/bits/rlevector.h incbwt/sasamples.h incbwt/misc/definitions.h \
+ incbwt/bits/bitbuffer.h incbwt/bits/deltavector.h incbwt/lcpsamples.h \
+ incbwt/bits/array.h incbwt/misc/utils.h incbwt/misc/definitions.h \
+ incbwt/misc/parameters.h RLCSABuilder.h TextCollectionBuilder.h \
+ TextCollection.h Tools.h TextStorage.h incbwt/bits/deltavector.h \
+ RLCSAWrapper.h incbwt/rlcsa.h
TextCollection.o: TextCollection.cpp TextCollection.h Tools.h FMIndex.h \
incbwt/bits/deltavector.h incbwt/bits/bitvector.h \
incbwt/bits/../misc/definitions.h incbwt/bits/bitbuffer.h BitRank.h \
BlockArray.h TextStorage.h ArrayDoc.h SWCSAWrapper.h \
- swcsa/utils/defValues.h swcsa/utils/valstring.h swcsa/interface.h
+ swcsa/utils/defValues.h swcsa/utils/valstring.h swcsa/interface.h \
+ RLCSAWrapper.h incbwt/rlcsa.h incbwt/bits/vectors.h \
+ incbwt/bits/deltavector.h incbwt/bits/rlevector.h incbwt/sasamples.h \
+ incbwt/misc/definitions.h incbwt/bits/bitbuffer.h \
+ incbwt/bits/deltavector.h incbwt/lcpsamples.h incbwt/bits/array.h \
+ incbwt/misc/utils.h incbwt/misc/definitions.h incbwt/misc/parameters.h
TextCollectionBuilder.o: TextCollectionBuilder.cpp \
TextCollectionBuilder.h TextCollection.h Tools.h TextStorage.h \
incbwt/bits/deltavector.h incbwt/bits/bitvector.h \
incbwt/bits/../misc/definitions.h incbwt/bits/bitbuffer.h \
FMIndexBuilder.h SWCSABuilder.h SWCSAWrapper.h swcsa/utils/defValues.h \
- swcsa/utils/valstring.h swcsa/interface.h
+ swcsa/utils/valstring.h swcsa/interface.h RLCSABuilder.h
TextStorage.o: TextStorage.cpp TextStorage.h TextCollection.h Tools.h \
incbwt/bits/deltavector.h incbwt/bits/bitvector.h \
incbwt/bits/../misc/definitions.h incbwt/bits/bitbuffer.h \
--- /dev/null
+Copyright (c) 2007-2009 Jouni Siren
+
+Permission is hereby granted, free of charge, to any person
+obtaining a copy of this software and associated documentation
+files (the "Software"), to deal in the Software without
+restriction, including without limitation the rights to use,
+copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following
+conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+OTHER DEALINGS IN THE SOFTWARE.
CC = g++
-CPPFLAGS = -Wall -O3
-OBJS = rlcsa.o rlcsa_builder.o sasamples.o bits/bitvector.o bits/deltavector.o bits/rlevector.o bits/vectors.o misc/parameters.o misc/utils.o qsufsort/qsufsort.o
+# Use 64-bit integers in a 64-bit environment.
+# SIZE_FLAGS = -DMASSIVE_DATA_RLCSA
-all: library
+# Parallelism is supported by either libstdc++ Parallel Mode or MCSTL.
+# PARALLEL_FLAGS = -DMULTITHREAD_SUPPORT -D_GLIBCXX_PARALLEL -fopenmp
+# MCSTL_ROOT = /fs-3/d/jltsiren/suds/mcstl
+# PARALLEL_FLAGS = -DMULTITHREAD_SUPPORT -I$(MCSTL_ROOT)/c++ -fopenmp
-library: $(OBJS)
- ar rc rlcsa.a $(OBJS)
+# Vectors using nibble codes instead of delta codes are faster, but they also
+# take up more space.
+# VECTOR_FLAGS = -DUSE_NIBBLE_VECTORS
+
+#CPPFLAGS = -Wall -g $(SIZE_FLAGS) $(PARALLEL_FLAGS) $(VECTOR_FLAGS)
+CPPFLAGS = -Wall -O3 $(SIZE_FLAGS) $(PARALLEL_FLAGS) $(VECTOR_FLAGS)
+OBJS = rlcsa.o rlcsa_builder.o sasamples.o lcpsamples.o bits/array.o bits/bitvector.o bits/deltavector.o bits/rlevector.o bits/nibblevector.o misc/parameters.o misc/utils.o qsufsort/qsufsort.o
+
+VPATH = bits:misc:utils
+
+
+default: rlcsa.a
+
+
+rlcsa.a: $(OBJS)
+ ar rcs rlcsa.a $(OBJS)
depend:
- g++ -MM *.cpp bits/*.cpp misc/*.cpp qsufsort/*.c > dependencies.mk
+ g++ -MM *.cpp bits/*.cpp misc/*.cpp qsufsort/*.c utils/*.cpp > dependencies.mk
+
+rlcsa_test: rlcsa_test.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o rlcsa_test rlcsa_test.o rlcsa.a
+
+lcp_test: lcp_test.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o lcp_test lcp_test.o rlcsa.a
+
+parallel_build: parallel_build.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o parallel_build parallel_build.o rlcsa.a
+
+build_rlcsa: build_rlcsa.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o build_rlcsa build_rlcsa.o rlcsa.a
+
+locate_test: locate_test.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o locate_test locate_test.o rlcsa.a
+
+display_test: display_test.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o display_test display_test.o rlcsa.a
+
+read_bwt: read_bwt.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o read_bwt read_bwt.o rlcsa.a
+
+extract_sequence: extract_sequence.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o extract_sequence extract_sequence.o rlcsa.a
+
+rlcsa_grep: rlcsa_grep.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o rlcsa_grep rlcsa_grep.o rlcsa.a
+
+build_plcp: build_plcp.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o build_plcp build_plcp.o rlcsa.a
+
+sample_lcp: sample_lcp.o rlcsa.a
+ $(CC) $(CPPFLAGS) -o sample_lcp sample_lcp.o rlcsa.a
+
+extract_text: extract_text.o
+ $(CC) $(CPPFLAGS) -o utils/extract_text extract_text.o
+
+convert_patterns: convert_patterns.o
+ $(CC) $(CPPFLAGS) -o utils/convert_patterns convert_patterns.o misc/utils.o
+
+split_text: split_text.o
+ $(CC) $(CPPFLAGS) -o utils/split_text split_text.o misc/utils.o
-update:
- cp ../rlcsa.* ../rlcsa_builder.* ../sasamples.* .
- cp ../bits/* bits/
- cp ../misc/* misc/
- cp ../qsufsort/* qsufsort/
+genpatterns: genpatterns.c
+ gcc -O3 -Wall -o utils/genpatterns utils/genpatterns.c
clean:
- rm -f rlcsa.a
- rm -f $(OBJS)
+ rm -f rlcsa.a rlcsa_test lcp_test parallel_build
+ rm -f build_rlcsa locate_test display_test read_bwt extract_sequence rlcsa_grep build_plcp sample_lcp
+ rm -f utils/extract_text utils/convert_patterns utils/split_text
+ rm -f *.o bits/*.o misc/*.o qsufsort/*.o utils/*.o
+package:
+ mkdir rlcsa
+ mkdir rlcsa/bits rlcsa/misc rlcsa/qsufsort rlcsa/utils
+ cp LICENSE Makefile README dependencies.mk *.cpp *.h rlcsa
+ cp bits/*.cpp bits/*.h rlcsa/bits
+ cp misc/*.cpp misc/*.h rlcsa/misc
+ cp qsufsort/*.c qsufsort/*.h rlcsa/qsufsort
+ cp utils/*.cpp utils/*.py rlcsa/utils
+ tar cfz rlcsa.tgz rlcsa
+ rm -r rlcsa/*
+ rmdir rlcsa
include dependencies.mk
--- /dev/null
+General Information
+===================
+
+This is an implementation of the Run-Length Compressed Suffix Array (RLCSA) [1,2] and its incremental construction [2]. The implementation includes experimental support for LCP information [3].
+
+Copyright 2007 - 2010, Jouni Siren, except when otherwise noted. See LICENSE for further information.
+
+This package contains an implementation of the algorithm presented in "Faster Suffix Sorting" by N. Jesper Larsson (jesper@cs.lth.se) and Kunihiko Sadakane (sada@is.s.u-tokyo.ac.jp). Copyright 1999, N. Jesper Larsson.
+
+
+Compiling
+---------
+
+The code should compile in both 32-bit and 64-bit environments. Uncomment SIZE_FLAGS in the makefile to use 64-bit integers in the 64-bit version.
+
+Parallelism is supported by libstdc++ Parallel Mode and by MCSTL. Uncomment either version of PARALLEL_FLAGS to compile the parallel version of the library, and set MCSTL_ROOT if necessary. GCC 4.2 or newer is required for the MCSTL version.
+
+Uncomment VECTOR_FLAGS to use a faster encoding for the run-length encoded bit vectors in .rlcsa.array. This increases the size somewhat.
+
+32-bit integers limit the size of the collection to less than 4 gigabytes. The size of individual input files is limited to less than 2 gigabytes in both 32-bit and 64-bit versions.
+
+Note that if 32-bit integers are used, then the bit-aligned arrays are limited to less than 512 megabytes (2^32 bits) in size. Hence if n is the collection size in characters and d is the sample rate, then (n / d) log ceil(n / d) must be less than 2^32. Otherwise the suffix array samples cannot be stored.
+
+
+Index Construction
+------------------
+
+The naming conventions for files are:
+
+ base_name - the sequences
+ base_name.rlcsa.array - most of the CSA
+ base_name.rlcsa.sa_samples - suffix array samples for locate and display
+ base_name.rlcsa.parameters - some of the index parameters
+ base_name.lcp_samples - sampled LCP array
+ base_name.plcp - run-length encoded PLCP array
+
+A typical parameter file looks like:
+
+ RLCSA_BLOCK_SIZE = 32
+ SAMPLE_RATE = 64
+ SUPPORT_DISPLAY = 1
+ SUPPORT_LOCATE = 1
+
+parallel_build is used to construct the index, as described in [2]. The program takes 2 to 4 parameters:
+
+ use '-n' as the first parameter to stop before merging the partial indexes
+ a list of input files (a text file, one file name per line)
+ base name of the output
+ number of threads to use (optional)
+
+The default parameters are 32 bytes for block size and 512 for sample rate. To modify these, one should create the parameter file before running the construction program. Each input file should be a concatenation of C-style '\0'-terminated strings. The files must be smaller than 2 GB each.
+
+build_rlcsa indexes a text file using a single thread. The program takes 2-4 parameters:
+
+ name of the input file
+ buffer size (in megabytes)
+ block size (optional)
+ sample rate (optional)
+
+Every line is inserted as a separate sequence. The lines must not contain character '\0'. The buffer size must be less than 2 gigabytes.
+
+
+Operations
+----------
+
+A number of operations have been implemented. The most important ones are the following:
+
+ pair_type count(const std::string& pattern)
+ Returns the suffix array range corresponding to the matches of the pattern. The range is reported as a closed interval.
+
+ usint* locate(pair_type range)
+ usint* locate(pair_type range, usint* data)
+ usint locate(usint index)
+ These return the suffix array values at given range/position. The user is responsible for the allocated data.
+
+ uchar* display(usint sequence)
+ uchar* display(usint sequence, pair_type range)
+ uchar* display(usint sequence, pair_type range, uchar* data)
+ These return a substring of the given sequence, as determined by the closed interval 'range'. The user is responsible for the allocated data.
+
+ uchar* display(usint position, usint len, usint context)
+ This is intended for displaying an occurrence of a pattern of length 'len' at position 'position' with 'context' extra characters on both sides.
+
+ uchar* readBWT()
+ uchar* readBWT(pair_type range)
+ Returns the BWT of the collection or a part of it. The user is responsible for the allocated string. Note that unlike the suffix array, the BWT includes all end markers.
+
+ pair_type getSequenceRange(usint number)
+ T^{number} -> [sp, ep], where T^{number} === T[sp, ep]
+
+ pair_type getSequenceRangeForPosition(usint value)
+ [sp, ep], where sp <= value <= ep
+
+ usint getSequenceForPosition(usint value)
+ i, where T^{i} === T[sp, ep] and sp <= value <= ep
+
+ usint* getSequenceForPosition(usint* value, usint length)
+ As above, but for multiple positions at once.
+
+ pair_type getRelativePosition(usint value)
+ (i, offset), where T^{i}[offset] === T[value]
+
+Locate and display can only be used when the corresponding parameter (SUPPORT_LOCATE or SUPPORT_DISPLAY) has value 1 and the suffix array samples have been created during the construction. If both of the parameters are missing or have value 0, the suffix array samples will not be loaded into memory.
+
+These operations are const and hence thread-safe.
+
+
+Construction Interface
+----------------------
+
+Class RLCSABuilder provides other possibilities for index construction. The constructor takes four parameters:
+
+ block size for the Psi vectors in bytes
+ suffix array sample rate
+ buffer size for construction in bytes
+ number of threads to use
+
+If suffix array sample rate is set to 0, the samples will not be created. The buffer size must be less than 2 gigabytes.
+
+Function insertSequence is called to insert a new sequence into the collection. The parameters are:
+
+ sequence as a char array
+ length of the sequence (not including the trailing 0, if present)
+ should we free the memory used by the sequence
+
+Function insertFromFile can be used to merge existing indexes into the collection. It takes the base name of the index as a parameter. The sequences and the index should both be available.
+
+Function getRLCSA is used to finish the construction and get the final index. After the call, the builder no longer contains the index. The caller is responsible for freeing the index.
+
+For example, the following inserts the sequences into the collection one at a time:
+
+ // Use a buffer of n megabytes.
+ RLCSABuilder builder(block_size, sample_rate, n * MEGABYTE, threads);
+
+ // For each sequence:
+ builder.insertSequence(sequence, length, false);
+
+ // If succesful, write the index to disk.
+ if(builder.isOk())
+ {
+ RLCSA* rlcsa = builder.getRLCSA();
+ rlcsa->writeTo(base_name);
+ delete rlcsa;
+ }
+
+
+Incremental construction for multiple sequences
+-----------------------------------------------
+
+When there are multiple sequences in one increment, character '\0' is assumed to represent the end of sequence marker. Hence the sequences themselves cannot contain character '\0'. This is always the case when using RLCSABuilder to build the partial indexes.
+
+If there is just one sequence in the increment, character '\0' is considered a normal character. This requires setting multiple_sequences = false in the RLCSA constructor. Note that RLCSABuilder cannot be used to merge these indexes, as it assumes character '\0' an end of sequence marker.
+
+
+LCP Support
+-----------
+
+The implementation includes experimental support for two representations of the LCP array: run-length encoded PLCP array and the sampled LCP array. sample_lcp and build_plcp can be used to build the representations. lcp_test was used in the experiments reported in [3].
+
+
+Other Programs
+--------------
+
+extract_sequence can be used to extract individual sequences from the index.
+
+rlcsa_grep is an example program offering limited grep-like functionality on indexes built with build_rlcsa.
+
+rlcsa_test and display_test were used in the experiments reported in [2]. Some small modifications have been required to use the current interface.
+
+The rest of the programs have not been used recently. They might no longer work correctly.
+
+
+Technical Information
+=====================
+
+A collection of sequences
+-------------------------
+
+The index contains a collection C of sequences T1, T2, ..., Tr. When constructing the index, each Ti is assumed to be followed by an implicit end of sequence marker $. The markers are assumed to be less than any character in the alphabet, and their mutual order is defined by the sequences they are following.
+
+The end of sequence markers are not included in the sequences or their suffix arrays. They are implicitly included in the values of Psi array as the first r values. If Psi(i) = j, then the pointer to the suffix starting at position SA[i] + 1 is stored at SA[Psi(i) - r]. These values are not stored in the index, however, making the Psi array effectively a collection of cycles instead of a single cycle.
+
+We use a bit vector E to mark the ending positions of each sequence. The indexed sequence is implicitly padded with empty characters so that each sequence starts at a position divisible by sample rate d. The starting position of sequence k > 0 is d * ((E.select(k - 1) / d) + 1).
+
+
+Suffix array samples
+--------------------
+
+For a given sample rate d, we store the positions of each sequence divisible by d. The sampled positions of suffix array are marked in a bit vector S, while the multipliers of d are stored in an array A in the same order. Another array B contains the inverse permutation of A.
+
+When locating, we can use S.valueAfter(i - r) to get (j, k = S.rank(j) - 1) for the first sampled j >= i - r in the suffix array order. If j == i - r, we can get the suffix array value as d * A[k]. If i < r, we have reached the implicit end of sequence marker. In this case, the suffix array value is E.select(i) + 1 (this value should only be used to derive SA values for earlier positions). Note that i is a value of Psi, not an index of the suffix array.
+
+When displaying text starting from position i, j = B[i / d] gives us the sample used as a starting point. The sample is located in position k = S.select(j) of suffix array corresponding to position k + r of Psi.
+
+
+Data formats
+------------
+
+.rlcsa.array
+ distribution of characters (CHARS * sizeof(usint) bytes)
+ RLEVector for each character appearing in the text
+ DeltaVector E
+ sample rate d (sizeof(usint) bytes)
+
+.rlcsa.sa_samples
+ DeltaVector S
+ array A (number_of_samples items of length(number_of_samples - 1) bits)
+
+Any bit vector
+ universe size (sizeof(usint) bytes)
+ item count (sizeof(usint) bytes)
+ number of blocks (sizeof(usint) bytes)
+ block size in words (sizeof(usint) bytes)
+ block data (number_of_blocks * block_size words)
+ sample array (2 * (number_of_blocks + 1) items of length(size) bits)
+
+
+LCP Information
+---------------
+
+The (P)LCP representations have been genralized to support multiple sequences. As the end markers are not included in the collection, the LCP values corresponding to the last characters of the sequences can be 1 or 0. The padding characters between the sequences are also assigned LCP values in the PLCP representation to ease its use. The sampled LCP array is used in a similar way as the SA samples in locate.
+
+Data formats:
+
+.lcp_samples
+ DeltaVector for the sampled positions
+ Array for the sampled values
+
+.plcp
+ DeltaVector
+
+Array
+ item count (sizeof(usint bytes))
+ number of blocks (sizeof(usint) bytes)
+ block size in words (sizeof(usint) bytes)
+ block data (number_of_blocks * block_size words)
+ file.write((char*)&(this->number_of_blocks), sizeof(this->number_of_blocks));
+ file.write((char*)&(this->block_size), sizeof(this->block_size));
+ file.write((char*)(this->array), this->block_size * this->number_of_blocks * sizeof(usint));
+ sample array (number_of_blocks + 1 items of length(items) bits)
+
+
+References
+==========
+
+[1] Veli Mäkinen, Gonzalo Navarro, Jouni Sirén, and Niko Välimäki: Storage and Retrieval of Highly Repetitive Sequence Collections.
+Journal of Computational Biology 17(3):281-308, 2010.
+
+[2] Jouni Sirén: Compressed Suffix Arrays for Massive Data.
+In SPIRE 2009, Springer LNCS 5721, pp. 63-74, Saariselkä, Finland, August 25-27, 2009.
+
+[3] Jouni Sirén: Sampled Longest Common Prefix Array.
+In CPM 2010, Springer LNCS 6129, pp. 227-237, New York, USA, June 21-23, 2010.
--- /dev/null
+#include <cstdlib>
+
+#include "array.h"
+
+
+namespace CSA
+{
+
+
+Array::Array(std::ifstream& file) :
+ delete_array(true), item_index(0)
+{
+ file.read((char*)&(this->items), sizeof(this->items));
+ file.read((char*)&(this->number_of_blocks), sizeof(this->number_of_blocks));
+ file.read((char*)&(this->block_size), sizeof(this->block_size));
+
+ this->array = new usint[this->block_size * this->number_of_blocks];
+ file.read((char*)(this->array), this->block_size * this->number_of_blocks * sizeof(usint));
+ this->buffer = new ReadBuffer(this->array, this->block_size);
+
+ this->integer_bits = length(this->items);
+ this->samples = new ReadBuffer(file, this->number_of_blocks + 1, this->integer_bits);
+
+ this->buildIndex();
+}
+
+Array::Array(ArrayEncoder& encoder) :
+ items(encoder.items),
+ delete_array(true),
+ block_size(encoder.block_size),
+ number_of_blocks(encoder.blocks),
+ item_index(0)
+{
+ usint* array_buffer = new usint[this->block_size * this->number_of_blocks];
+ this->array = array_buffer;
+ this->buffer = new ReadBuffer(this->array, this->block_size);
+
+ this->integer_bits = length(this->items);
+ WriteBuffer sample_buffer(this->number_of_blocks + 1, this->integer_bits);
+
+ // Copy & linearize the array.
+ usint pos = 0;
+ for(std::list<usint*>::iterator iter = encoder.array_blocks.begin(); iter != encoder.array_blocks.end(); iter++)
+ {
+ memcpy(array_buffer + pos, *iter, encoder.superblock_bytes);
+ pos += encoder.block_size * encoder.blocks_in_superblock;
+ }
+ memcpy(array_buffer + pos, encoder.array, encoder.current_blocks * encoder.block_size * sizeof(usint));
+
+ // Compress the samples.
+ for(std::list<usint*>::iterator iter = encoder.sample_blocks.begin(); iter != encoder.sample_blocks.end(); iter++)
+ {
+ usint* buf = *iter;
+ for(usint i = 0; i < encoder.samples_in_superblock; i++)
+ {
+ sample_buffer.writeItem(buf[i]);
+ }
+ }
+ for(usint i = 0; i < encoder.current_samples; i++)
+ {
+ sample_buffer.writeItem(encoder.samples[i]);
+ }
+ sample_buffer.writeItem(this->items);
+
+ this->samples = sample_buffer.getReadBuffer();
+
+ this->buildIndex();
+}
+
+Array::~Array()
+{
+ if(this->delete_array) { delete[] this->array; }
+ delete this->buffer;
+ delete this->samples;
+ delete this->item_index;
+}
+
+//--------------------------------------------------------------------------
+
+void
+Array::writeTo(std::ofstream& file) const
+{
+ file.write((char*)&(this->items), sizeof(this->items));
+ file.write((char*)&(this->number_of_blocks), sizeof(this->number_of_blocks));
+ file.write((char*)&(this->block_size), sizeof(this->block_size));
+ file.write((char*)(this->array), this->block_size * this->number_of_blocks * sizeof(usint));
+ this->samples->writeBuffer(file);
+}
+
+usint
+Array::reportSize() const
+{
+ usint bytes = sizeof(*this);
+ bytes += this->buffer->reportSize();
+ bytes += this->block_size * this->number_of_blocks * sizeof(usint);
+ bytes += this->samples->reportSize();
+ bytes += this->item_index->reportSize();
+ return bytes;
+}
+
+//--------------------------------------------------------------------------
+
+void
+Array::buildIndex()
+{
+ delete this->item_index;
+
+ usint index_samples = (this->number_of_blocks + Array::INDEX_RATE - 1) / Array::INDEX_RATE;
+ this->index_rate = (this->items + index_samples - 1) / index_samples;
+ index_samples = (this->items + this->index_rate - 1) / this->index_rate + 1;
+ WriteBuffer index_buffer(index_samples, length(this->number_of_blocks - 1));
+
+ usint current = 0, pointer = 0;
+ this->samples->goToItem(1);
+ while(this->samples->hasNextItem())
+ {
+ usint limit = this->samples->readItem();
+ while(current < limit)
+ {
+ index_buffer.writeItem(pointer);
+ current += this->index_rate;
+ }
+ pointer++;
+ }
+ index_buffer.writeItem(this->number_of_blocks - 1);
+
+ this->item_index = index_buffer.getReadBuffer();
+}
+
+//--------------------------------------------------------------------------
+
+Array::Iterator::Iterator(const Array& par) :
+ parent(par),
+ buffer(*(par.buffer)),
+ samples(*(par.samples))
+{
+}
+
+Array::Iterator::~Iterator()
+{
+}
+
+usint
+Array::Iterator::getItem(usint index)
+{
+ if(index >= this->parent.items) { return 0; }
+ this->getSample(this->sampleForIndex(index));
+
+ usint value = 0;
+ while(this->sample + this->cur <= index)
+ {
+ value = this->buffer.readDeltaCode();
+ this->cur++;
+ }
+
+ return value;
+}
+
+usint
+Array::Iterator::nextItem()
+{
+ if(this->cur >= this->block_items)
+ {
+ this->getSample(this->block + 1);
+ }
+
+ this->cur++;
+ return this->buffer.readDeltaCode();
+}
+
+usint
+Array::Iterator::sampleForIndex(usint index)
+{
+ usint low = this->parent.item_index->readItemConst(index / this->parent.index_rate);
+ usint high = this->parent.number_of_blocks - 1;
+
+ this->samples.goToItem(low + 1);
+ for(; low < high; low++)
+ {
+ if(this->samples.readItem() > index) { return low; }
+ }
+
+ return low;
+}
+
+//--------------------------------------------------------------------------
+
+ArrayEncoder::ArrayEncoder(usint block_bytes, usint superblock_size) :
+ items(0), blocks(1),
+ block_size(BYTES_TO_WORDS(block_bytes)),
+ superblock_bytes(superblock_size)
+{
+ this->array = new usint[this->superblock_bytes / sizeof(usint)];
+ memset(this->array, 0, this->superblock_bytes);
+ this->blocks_in_superblock = this->superblock_bytes / (sizeof(usint) * this->block_size);
+ this->current_blocks = 1;
+
+ this->samples = new usint[this->superblock_bytes / sizeof(usint)];
+ this->samples_in_superblock = this->superblock_bytes / sizeof(usint);
+ this->current_samples = 1;
+ this->samples[0] = 0;
+
+ this->buffer = new WriteBuffer(this->array, this->block_size);
+}
+
+ArrayEncoder::~ArrayEncoder()
+{
+ delete[] this->array;
+
+ delete this->buffer;
+ for(std::list<usint*>::iterator iter = this->array_blocks.begin(); iter != this->array_blocks.end(); iter++)
+ {
+ delete[] *iter;
+ }
+
+ delete[] this->samples;
+ for(std::list<usint*>::iterator iter = this->sample_blocks.begin(); iter != this->sample_blocks.end(); iter++)
+ {
+ delete[] *iter;
+ }
+}
+
+void
+ArrayEncoder::addItem(usint value)
+{
+ this->items++;
+ if(this->buffer->writeDeltaCode(value)) { return; }
+
+ // Didn't fit into the block. A new block & sample required.
+ this->blocks++;
+ this->current_blocks++;
+ this->current_samples++;
+
+ // Do we need a new superblock for the block?
+ if(this->current_blocks > this->blocks_in_superblock)
+ {
+ this->array_blocks.push_back(this->array);
+ this->array = new usint[this->superblock_bytes / sizeof(usint)];
+ memset(this->array, 0, this->superblock_bytes);
+ this->current_blocks = 1;
+ }
+ this->buffer->moveBuffer(this->array + (this->block_size * (this->current_blocks - 1)));
+
+ // Do we need a new superblock for the sample?
+ if(this->current_samples > this->samples_in_superblock)
+ {
+ this->sample_blocks.push_back(this->samples);
+ this->samples = new usint[this->superblock_bytes / sizeof(usint)];
+ this->current_samples = 1;
+ }
+ this->samples[this->current_samples - 1] = this->items - 1;
+
+ // Finally, write the item.
+ this->buffer->writeDeltaCode(value);
+}
+
+
+} // namespace CSA
--- /dev/null
+#ifndef ARRAY_H
+#define ARRAY_H
+
+#include <fstream>
+#include <list>
+
+#include "../misc/definitions.h"
+#include "bitbuffer.h"
+
+
+namespace CSA
+{
+
+/*
+ Delta-coded array for storing positive integers.
+ Do not try to encode value 0!
+*/
+
+
+class ArrayEncoder
+{
+ public:
+ const static usint SUPERBLOCK_SIZE = MEGABYTE;
+
+ // We assume superblock size is divisible by block and sample size.
+ ArrayEncoder(usint block_bytes, usint superblock_size = SUPERBLOCK_SIZE);
+ ~ArrayEncoder();
+
+ void addItem(usint value);
+
+ void addNewBlock();
+
+ usint items, blocks;
+ usint block_size, superblock_bytes;
+
+ WriteBuffer* buffer;
+
+ std::list<usint*> array_blocks;
+ usint* array;
+ usint blocks_in_superblock, current_blocks;
+
+ std::list<usint*> sample_blocks;
+ usint* samples;
+ usint samples_in_superblock, current_samples;
+
+ protected:
+
+ // These are not allowed.
+ ArrayEncoder();
+ ArrayEncoder(const ArrayEncoder&);
+ ArrayEncoder& operator = (const ArrayEncoder&);
+};
+
+
+class Array
+{
+ public:
+ const static usint INDEX_RATE = 5;
+
+ Array(std::ifstream& file);
+ Array(ArrayEncoder& encoder);
+ ~Array();
+
+//--------------------------------------------------------------------------
+
+ void writeTo(std::ofstream& file) const;
+
+ inline usint getSize() const { return this->items; }
+ inline usint getBlockSize() const { return this->block_size; }
+
+ usint reportSize() const;
+
+//--------------------------------------------------------------------------
+
+ class Iterator
+ {
+ public:
+ Iterator(const Array& par);
+ ~Iterator();
+
+ // Returns 0 if the index is invalid.
+ usint getItem(usint index);
+ usint nextItem();
+
+ inline bool hasNext() const
+ {
+ return (this->sample + this->cur < this->parent.getSize());
+ }
+
+ private:
+ const Array& parent;
+
+ usint block, sample, cur, block_items;
+
+ ReadBuffer buffer, samples;
+
+ /*
+ This function returns the sample corresponding to the block the given
+ index might be found in. Parameter is assumed to be valid.
+ */
+ usint sampleForIndex(usint index);
+
+ inline void getSample(usint sample_number)
+ {
+ this->block = sample_number;
+ this->sample = this->samples.readItem(sample_number);
+ this->cur = 0;
+ this->block_items = this->samples.readItem() - this->sample;
+ this->buffer.moveBuffer(this->parent.array + (this->block * this->parent.block_size));
+ }
+
+ // These are not allowed.
+ Iterator();
+ Iterator(const Iterator&);
+ Iterator& operator = (const Iterator&);
+ };
+
+//--------------------------------------------------------------------------
+
+ protected:
+ usint items;
+
+ ReadBuffer* buffer;
+ const usint* array;
+ bool delete_array; // Do we own the array?
+ usint block_size;
+ usint number_of_blocks;
+
+ /*
+ Each sample tells how many items are contained in the previous blocks.
+ */
+ ReadBuffer* samples;
+ usint integer_bits;
+
+ ReadBuffer* item_index;
+ usint index_rate;
+
+ /*
+ Build a higher level index for faster queries.
+ The index consists of about (number of samples) / INDEX_RATE pointers.
+ The array cannot be used without the index.
+ */
+ void buildIndex();
+
+ // These are not allowed.
+ Array();
+ Array(const Array&);
+ Array& operator = (const Array&);
+};
+
+
+} // namespace CSA
+
+
+#endif // ARRAY_H
#include <algorithm>
#include <cstdlib>
+#include <cstring>
#include <fstream>
#include <iostream>
-#include <cstring> // defines std::memset, added by Kim
#include "../misc/definitions.h"
{
-template<class Data>
-class GenericBitBuffer
+class ReadBuffer
{
public:
- GenericBitBuffer(usint words);
- GenericBitBuffer(usint _items, usint item_size);
- GenericBitBuffer(std::ifstream& file, usint words);
- GenericBitBuffer(std::ifstream& file, usint _items, usint item_size);
- GenericBitBuffer(std::FILE* file, usint _items, usint item_size);
- GenericBitBuffer(Data* buffer, usint words);
- GenericBitBuffer(Data* buffer, usint _items, usint item_size);
- ~GenericBitBuffer();
-
- void writeBuffer(std::ofstream& file, bool erase = true);
- void writeBuffer(std::FILE * file, bool erase = true);
- void readBuffer(std::ifstream& file, usint words, bool erase = true);
- void setBuffer(Data* buffer, usint words);
-
- // We assume we are already using a buffer not owned by this.
- inline void moveBuffer(Data* buffer)
+ ReadBuffer(std::ifstream& file, usint words) :
+ size(words),
+ item_bits(1),
+ items(0),
+ free_buffer(true)
{
+ usint* buffer = new usint[this->size];
+ memset(buffer, 0, this->size * sizeof(usint));
+ file.read((char*)buffer, this->size * sizeof(usint));
this->data = buffer;
- this->pos = 0;
- this->bits = DATA_BITS;
- this->current = 0;
+ this->reset();
}
-//--------------------------------------------------------------------------
+ ReadBuffer(std::ifstream& file, usint _items, usint item_size) :
+ item_bits(item_size),
+ items(_items),
+ free_buffer(true)
+ {
+ this->size = bitsToWords(this->items * this->item_bits);
+ usint* buffer = new usint[this->size];
+ memset(buffer, 0, this->size * sizeof(usint));
+ file.read((char*)buffer, this->size * sizeof(usint));
+ this->data = buffer;
+ this->reset();
+ }
- inline void reset(bool erase)
+ ReadBuffer(std::FILE* file, usint _items, usint item_size) :
+ item_bits(item_size),
+ items(_items),
+ free_buffer(true)
{
- this->pos = 0;
- this->bits = DATA_BITS;
- this->current = 0;
- if(erase)
- {
- memset(this->data, 0, this->size * sizeof(Data));
- }
+ this->size = bitsToWords(this->items * this->item_bits);
+ usint* buffer = new usint[this->size];
+ memset(buffer, 0, this->size * sizeof(usint));
+ fread(buffer, sizeof(usint), this->size, file);
+ this->data = buffer;
+ this->reset();
}
- inline void skipBits(usint count)
+ // This version does not delete the data when deleted.
+ ReadBuffer(const usint* buffer, usint words) :
+ size(words),
+ item_bits(1),
+ items(0),
+ free_buffer(false)
{
- if(count < this->bits)
- {
- this->bits -= count;
- return;
- }
+ this->data = buffer;
+ this->reset();
+ }
- count -= this->bits;
- this->pos += 1 + count / DATA_BITS;
- this->bits = DATA_BITS - count % DATA_BITS;
+ // This version does not delete the data when deleted.
+ ReadBuffer(const usint* buffer, usint _items, usint item_size) :
+ item_bits(item_size),
+ items(_items),
+ free_buffer(false)
+ {
+ this->size = bitsToWords(this->items * this->item_bits);
+ this->data = buffer;
+ this->reset();
+ }
+
+ // This version does not delete the data when deleted.
+ ReadBuffer(const ReadBuffer& original) :
+ data(original.data),
+ size(original.size),
+ item_bits(original.item_bits),
+ items(original.items),
+ free_buffer(false)
+ {
+ this->reset();
}
- inline void rewind(usint count)
+ ~ReadBuffer()
{
- this->bits += count;
- if(this->bits > DATA_BITS)
+ if(this->free_buffer)
{
- usint back = (this->bits - 1) / DATA_BITS;
- this->pos -= back;
- this->bits -= back * DATA_BITS;
+ delete[] this->data;
}
}
//--------------------------------------------------------------------------
- inline usint bitsLeft()
+ void claimData()
{
- return this->bits + DATA_BITS * (this->size - this->pos - 1);
+ this->free_buffer = true;
}
- inline void writeBits(usint value, usint count)
+ void writeBuffer(std::ofstream& file) const
+ {
+ file.write((const char*)this->data, this->size * sizeof(usint));
+ }
+ void writeBuffer(std::FILE* file) const
{
- while(count >= this->bits)
+ fwrite(this->data, sizeof(usint), this->size, file);
+ }
+
+ // The buffer will no longer own the data.
+ void moveBuffer(const usint* buffer)
+ {
+ if(this->free_buffer)
{
- count -= this->bits;
- this->data[this->pos] |= GET(LOWER(value, count), this->bits);
- this->pos++; this->bits = DATA_BITS;
+ delete[] this->data;
}
- if(count > 0)
+ this->free_buffer = false;
+
+ this->data = buffer;
+ this->reset();
+ }
+
+ usint reportSize() const
+ {
+ usint bytes = sizeof(*this);
+ if(this->free_buffer) { bytes += this->size * sizeof(usint); }
+ return bytes;
+ }
+
+//--------------------------------------------------------------------------
+
+ inline void reset()
+ {
+ this->pos = 0;
+ this->bits = WORD_BITS;
+ this->current = 0;
+ }
+
+ inline void skipBits(usint count)
+ {
+ if(count < this->bits)
{
this->bits -= count;
- this->data[this->pos] |= HIGHER(GET(value, count), this->bits);
+ return;
}
+
+ count -= this->bits;
+ this->pos += 1 + count / WORD_BITS;
+ this->bits = WORD_BITS - count % WORD_BITS;
+ }
+
+//--------------------------------------------------------------------------
+
+ inline usint bitsLeft() const
+ {
+ return this->bits + WORD_BITS * (this->size - this->pos - 1);
}
// Returns nonzero if bit is 1
this->bits--;
usint bit = this->data[this->pos] & ((usint)1 << this->bits);
- if(this->bits == 0) { this->pos++; this->bits = DATA_BITS; }
+ if(this->bits == 0) { this->pos++; this->bits = WORD_BITS; }
return bit;
}
{
usint value = 0;
- while(count >= this->bits)
+ if(count >= this->bits)
{
count -= this->bits;
value |= HIGHER(GET(this->data[this->pos], this->bits), count);
- this->pos++; this->bits = DATA_BITS;
+ this->pos++; this->bits = WORD_BITS;
}
if(count > 0)
{
//--------------------------------------------------------------------------
/*
- These operations work on fixed-size items.
+ These operations work on 4-bit nibbles.
+ Do not use these with the bit-level operations.
*/
- inline usint getItemSize()
+ inline usint readNibble()
{
- return this->item_bits;
+ this->bits -= 4;
+ usint value = GET(LOWER(this->data[this->pos], this->bits), 4);
+
+ if(this->bits == 0) { this->pos++; this->bits = WORD_BITS; }
+
+ return value;
}
-/*
- inline void setItemSize(usint new_size)
+ // Nibble code for positive integers.
+ inline usint readNibbleCode()
{
- if(new_size == 0) { return; }
- this->item_bits = new_size;
+ usint temp, value = 0, shift = 0;
+ do
+ {
+ temp = this->readNibble();
+ value |= (temp & 0x7) << shift;
+ shift += 3;
+ }
+ while((temp & 0x8) == 0);
+
+ return value + 1;
+ }
+
+//--------------------------------------------------------------------------
+
+ /*
+ These operations work on fixed-size items. No sanity checks are made
+ for parameter values.
+ */
+
+ inline usint getItemSize() const
+ {
+ return this->item_bits;
}
-*/
inline void goToItem(usint item)
{
usint b = item * this->item_bits;
- this->pos = b / DATA_BITS;
- this->bits = DATA_BITS - b % DATA_BITS;
+ this->pos = b / WORD_BITS;
+ this->bits = WORD_BITS - b % WORD_BITS;
this->current = item;
}
return this->readItem(0);
}
- inline bool hasNextItem()
+ inline usint readItemConst(usint item) const
{
- return (this->current < this->items);
+ usint b = item * this->item_bits;
+ usint p = b / WORD_BITS;
+ b = WORD_BITS - b % WORD_BITS;
+
+ usint c = this->item_bits;
+ usint value = 0;
+
+ if(c >= b)
+ {
+ c -= b;
+ value |= HIGHER(GET(this->data[p], b), c);
+ p++; b = WORD_BITS;
+ }
+ if(c > 0)
+ {
+ b -= c;
+ value |= GET(LOWER(this->data[p], b), c);
+ }
+
+ return value;
}
- inline void writeItem(usint item)
+ inline bool hasNextItem() const
{
- this->writeBits(item, this->item_bits);
- this->current++;
+ return (this->current < this->items);
}
inline void skipItem()
Delta coding for positive integers
*/
- inline bool canDeltaCode(usint value)
+ inline usint readDeltaCode()
{
- return this->deltaCodeLength(value) <= this->bitsLeft();
+ usint len = 0;
+ while(this->readBit() == 0) { len++; }
+
+ usint temp = (((usint)1 << len) | this->readBits(len)) - 1;
+ temp = ((usint)1 << temp) | this->readBits(temp);
+ return temp;
}
- inline usint deltaCodeLength(usint value)
+//--------------------------------------------------------------------------
+
+ private:
+ const usint* data;
+ usint size, item_bits, items;
+ bool free_buffer;
+
+ // Iterator data
+ usint pos, bits, current;
+
+ inline static usint bitsToWords(usint _bits) { return (_bits + WORD_BITS - 1) / WORD_BITS; }
+
+ // These are not allowed.
+ ReadBuffer();
+ ReadBuffer& operator = (const ReadBuffer&);
+};
+
+
+//--------------------------------------------------------------------------
+
+
+class WriteBuffer
+{
+ public:
+ WriteBuffer(usint words) :
+ size(words),
+ item_bits(1),
+ items(0),
+ free_buffer(true)
{
- usint len = length(value);
- usint llen = length(len);
- return (len + llen + llen - 2);
+ this->data = new usint[words];
+ memset(this->data, 0, this->size * sizeof(usint));
+ this->reset();
}
- // This version returns false if there is no space left for the encoding.
- inline bool writeDeltaCode(usint value)
+ WriteBuffer(usint _items, usint item_size) :
+ item_bits(item_size),
+ items(_items),
+ free_buffer(true)
{
- usint len = length(value);
- usint llen = length(len);
-
- if(len + llen + llen - 2 > this->bitsLeft()) { return false; }
+ this->size = bitsToWords(this->items * this->item_bits);
+ this->data = new usint[this->size];
+ memset(this->data, 0, this->size * sizeof(usint));
+ this->reset();
+ }
- // this->writeBits(0, llen - 1); // Now included in the next writeBits()
- this->writeBits(len, llen + llen - 1);
- this->writeBits(value, len - 1);
- return true;
+ // This version does not delete the data when deleted.
+ WriteBuffer(usint* buffer, usint words) :
+ size(words),
+ item_bits(1),
+ items(0),
+ free_buffer(false)
+ {
+ this->data = buffer;
+ this->reset();
}
- // This version assumes the code fits into the buffer.
- inline void writeDeltaCodeDirect(usint value)
+ // This version does not delete the data when deleted.
+ WriteBuffer(usint* buffer, usint _items, usint item_size) :
+ item_bits(item_size),
+ items(_items),
+ free_buffer(false)
{
- usint len = length(value);
- usint llen = length(len);
+ this->size = bitsToWords(this->items * this->item_bits);
+ this->data = buffer;
+ this->reset();
+ }
- // this->writeBits(0, llen - 1); // Now included in the next writeBits()
- this->writeBits(len, llen + llen - 1);
- this->writeBits(value, len - 1);
+ ~WriteBuffer()
+ {
+ if(this->free_buffer)
+ {
+ delete[] this->data;
+ }
}
- // We assume the code fits into usint:
- // 32-bit: value < 2^24
- // 64-bit: value < 2^54
- inline void writeDeltaCodeFast(usint value)
+//--------------------------------------------------------------------------
+
+ // This transfers the ownership of the data to the read buffer.
+ ReadBuffer* getReadBuffer()
{
- usint len = length(value);
+ ReadBuffer* buffer;
+ if(this->items > 0)
+ {
+ buffer = new ReadBuffer(this->data, this->items, this->item_bits);
+ }
+ else
+ {
+ buffer = new ReadBuffer(this->data, this->size);
+ }
- value ^= ((usint)1 << (len - 1));
- this->writeBits((len << (len - 1)) | value, len + 2 * length(len) - 2);
+ if(this->free_buffer)
+ {
+ buffer->claimData();
+ this->free_buffer = false;
+ }
+
+ return buffer;
}
- inline usint readDeltaCode()
+ void writeBuffer(std::ofstream& file) const
{
- usint len = 0;
- while(this->readBit() == 0) { len++; }
+ file.write((char*)this->data, this->size * sizeof(usint));
+ }
+ void writeBuffer(std::FILE *file) const
+ {
+ fwrite(this->data, sizeof(usint), this->size, file);
+ }
- usint temp = (((usint)1 << len) | this->readBits(len)) - 1;
- temp = ((usint)1 << temp) | this->readBits(temp);
- return temp;
+ // The buffer will no longer own the data.
+ void moveBuffer(usint* buffer)
+ {
+ if(this->free_buffer)
+ {
+ delete[] this->data;
+ }
+ this->free_buffer = false;
+
+ this->data = buffer;
+ this->reset();
+ }
+
+ usint reportSize() const
+ {
+ usint bytes = sizeof(*this);
+ if(this->free_buffer) { bytes += this->size * sizeof(usint); }
+ return bytes;
}
//--------------------------------------------------------------------------
- /*
- Gamma coding for positive integers
- */
+ inline void reset()
+ {
+ this->pos = 0;
+ this->bits = WORD_BITS;
+ this->current = 0;
+ }
- inline bool canGammaCode(usint value)
+ inline void skipBits(usint count)
{
- return this->gammaCodeLength(value) <= this->bitsLeft();
+ if(count < this->bits)
+ {
+ this->bits -= count;
+ return;
+ }
+
+ count -= this->bits;
+ this->pos += 1 + count / WORD_BITS;
+ this->bits = WORD_BITS - count % WORD_BITS;
}
- inline usint gammaCodeLength(usint value)
+//--------------------------------------------------------------------------
+
+ inline usint bitsLeft() const
{
- return 2 * length(value) - 1;
+ return this->bits + WORD_BITS * (this->size - this->pos - 1);
}
- // This version returns false if there is no space left for the encoding.
- inline bool writeGammaCode(usint value)
+ inline void writeBits(usint value, usint count)
{
- usint len = length(value);
+ if(count >= this->bits)
+ {
+ count -= this->bits;
+ this->data[this->pos] |= GET(LOWER(value, count), this->bits);
+ this->pos++; this->bits = WORD_BITS;
+ }
+ if(count > 0)
+ {
+ this->bits -= count;
+ this->data[this->pos] |= HIGHER(GET(value, count), this->bits);
+ }
+ }
- if(len > this->bitsLeft()) { return false; }
+//--------------------------------------------------------------------------
- this->writeBits(0, len - 1);
- this->writeBits(value, len);
- return true;
+ /*
+ These operations work on fixed-size items.
+ */
+
+ inline usint getItemSize() const
+ {
+ return this->item_bits;
}
- // This version assumes the code fits into the buffer.
- inline void writeGammaCodeDirect(usint value)
+ inline void goToItem(usint item)
{
- usint len = length(value);
+ usint b = item * this->item_bits;
+ this->pos = b / WORD_BITS;
+ this->bits = WORD_BITS - b % WORD_BITS;
+ this->current = item;
+ }
- this->writeBits(0, len - 1);
- this->writeBits(value, len);
+ inline bool hasNextItem() const
+ {
+ return (this->current < this->items);
}
- // We assume the code fits into usint:
- // 32-bit: value < 2^16
- // 64-bit: value < 2^32
- inline void writeGammaCodeFast(usint value)
+ inline void writeItem(usint item)
{
- this->writeBits(value, this->gammaCodeLength(value));
+ this->writeBits(item, this->item_bits);
+ this->current++;
}
- inline usint readGammaCode()
+ inline void skipItem()
{
- usint len = 1;
- while(this->readBit() == 0) { len++; }
- return ((usint)1 << len) | this->readBits(len);
+ this->skipBits(this->item_bits);
+ this->current++;
}
//--------------------------------------------------------------------------
- usint reportSize();
+ /*
+ Nibble coding for positive integers.
+ */
-//--------------------------------------------------------------------------
+ inline usint nibbleCodeLength(usint value) const
+ {
+ usint b = 0;
+ value--;
- private:
- Data* data;
- usint size, pos, bits;
- usint item_bits, items, current;
- bool free_buffer;
+ do
+ {
+ b += 4;
+ value >>= 3;
+ }
+ while(value > 0);
- const static usint DATA_BITS = sizeof(Data) * CHAR_BIT;
+ return b;
+ }
- inline static usint bitsToData(usint _bits) { return (_bits + DATA_BITS - 1) / DATA_BITS; }
-};
+ // Something breaks very badly if value > 15.
+ inline void writeNibble(usint value)
+ {
+ this->bits -= 4;
+ this->data[this->pos] |= HIGHER(value, this->bits);
+ if(this->bits == 0) { this->pos++; this->bits = WORD_BITS; }
+ }
+ // It is assumed that there is enough space for the code.
+ inline void writeNibbleCode(usint value)
+ {
+ value--;
+ while(value > 0x7)
+ {
+ this->writeNibble(value & 0x7);
+ value >>= 3;
+ }
+ this->writeNibble(value | 0x8);
+ }
-typedef GenericBitBuffer<uchar> BitBuffer;
-typedef GenericBitBuffer<usint> FastBitBuffer;
+//--------------------------------------------------------------------------
+ /*
+ Delta coding for positive integers
+ */
-//--------------------------------------------------------------------------
+ inline bool canDeltaCode(usint value) const
+ {
+ return this->deltaCodeLength(value) <= this->bitsLeft();
+ }
+ inline usint deltaCodeLength(usint value) const
+ {
+ usint len = length(value);
+ usint llen = length(len);
+ return (len + llen + llen - 2);
+ }
-template<class Data>
-GenericBitBuffer<Data>::GenericBitBuffer(usint words) :
- size(words),
- pos(0),
- bits(DATA_BITS),
- item_bits(1),
- items(0),
- current(0),
- free_buffer(true)
-{
- this->data = new Data[words];
- memset(this->data, 0, words * sizeof(Data));
-}
-
-template<class Data>
-GenericBitBuffer<Data>::GenericBitBuffer(usint _items, usint item_size) :
- pos(0),
- bits(DATA_BITS),
- item_bits(item_size),
- items(_items),
- current(0),
- free_buffer(true)
-{
- this->size = bitsToData(items * item_size);
- this->data = new Data[this->size];
- memset(this->data, 0, this->size * sizeof(Data));
-}
-
-template<class Data>
-GenericBitBuffer<Data>::GenericBitBuffer(std::ifstream& file, usint words) :
- size(words),
- pos(0),
- bits(DATA_BITS),
- item_bits(1),
- items(0),
- current(0),
- free_buffer(true)
-{
- this->data = new Data[words];
- memset(this->data, 0, words * sizeof(Data));
- file.read((char*)this->data, words * sizeof(Data));
-}
-
-template<class Data>
-GenericBitBuffer<Data>::GenericBitBuffer(std::ifstream& file, usint _items, usint item_size) :
- size(BITS_TO_WORDS(_items * item_size)),
- pos(0),
- bits(DATA_BITS),
- item_bits(item_size),
- items(_items),
- current(0),
- free_buffer(true)
-{
- this->data = new Data[this->size];
- memset(this->data, 0, this->size * sizeof(Data));
- file.read((char*)this->data, this->size * sizeof(Data));
-}
-
-template<class Data>
-GenericBitBuffer<Data>::GenericBitBuffer(std::FILE* file, usint _items, usint item_size) :
- size(BITS_TO_WORDS(_items * item_size)),
- pos(0),
- bits(DATA_BITS),
- item_bits(item_size),
- items(_items),
- current(0),
- free_buffer(true)
-{
- this->data = new Data[this->size];
- memset(this->data, 0, this->size * sizeof(Data));
- std::fread(this->data, sizeof(Data), this->size, file);
-}
-
-template<class Data>
-GenericBitBuffer<Data>::GenericBitBuffer(Data* buffer, usint words) :
- size(words),
- pos(0),
- bits(DATA_BITS),
- item_bits(1),
- items(0),
- current(0),
- free_buffer(false)
-{
- this->data = buffer;
-}
-
-template<class Data>
-GenericBitBuffer<Data>::GenericBitBuffer(Data* buffer, usint _items, usint item_size) :
- size(BITS_TO_WORDS(_items * item_size)),
- pos(0),
- bits(DATA_BITS),
- item_bits(item_size),
- items(_items),
- current(0),
- free_buffer(false)
-{
- this->data = buffer;
-}
+ // This version returns false if there is no space left for the encoding.
+ inline bool writeDeltaCode(usint value)
+ {
+ usint len = length(value);
+ usint llen = length(len);
-template<class Data>
-GenericBitBuffer<Data>::~GenericBitBuffer()
-{
- if(this->free_buffer)
- {
- delete[] this->data;
- }
-}
+ if(len + llen + llen - 2 > this->bitsLeft()) { return false; }
-//--------------------------------------------------------------------------
+ // this->writeBits(0, llen - 1); // Now included in the next writeBits()
+ this->writeBits(len, llen + llen - 1);
+ this->writeBits(value, len - 1);
+ return true;
+ }
-template<class Data>
-void
-GenericBitBuffer<Data>::writeBuffer(std::ofstream& file, bool erase)
-{
- file.write((char*)this->data, this->size * sizeof(Data));
- this->reset(erase);
-}
+ // This version assumes the code fits into the buffer.
+ inline void writeDeltaCodeDirect(usint value)
+ {
+ usint len = length(value);
+ usint llen = length(len);
-template<class Data>
-void
-GenericBitBuffer<Data>::writeBuffer(std::FILE* file, bool erase)
-{
- std::fwrite(this->data, sizeof(Data), this->size, file);
- this->reset(erase);
-}
+ // this->writeBits(0, llen - 1); // Now included in the next writeBits()
+ this->writeBits(len, llen + llen - 1);
+ this->writeBits(value, len - 1);
+ }
-template<class Data>
-void
-GenericBitBuffer<Data>::readBuffer(std::ifstream& file, usint words, bool erase)
-{
- if(words > this->size || !(this->free_buffer))
- {
- if(this->free_buffer)
+ // We assume the code fits into usint:
+ // 32-bit: value < 2^24
+ // 64-bit: value < 2^54
+ inline void writeDeltaCodeFast(usint value)
{
- delete[] this->data;
+ usint len = length(value);
+
+ value ^= ((usint)1 << (len - 1));
+ this->writeBits((len << (len - 1)) | value, len + 2 * length(len) - 2);
}
- this->size = words;
- this->data = new Data[words];
- this->free_buffer = true;
- }
- this->reset(erase);
- file.read((char*)this->data, words * sizeof(Data));
-}
+//--------------------------------------------------------------------------
-template<class Data>
-void
-GenericBitBuffer<Data>::setBuffer(Data* buffer, usint words)
-{
- if(this->free_buffer)
- {
- delete[] this->data;
- this->free_buffer = false;
- }
+ private:
+ usint* data;
+ usint size, item_bits, items;
+ bool free_buffer;
- this->data = buffer;
- this->size = words;
- this->reset(false);
-}
+ // Iterator data
+ usint pos, bits, current;
-//--------------------------------------------------------------------------
+ inline static usint bitsToWords(usint _bits) { return (_bits + WORD_BITS - 1) / WORD_BITS; }
+
+ // These are not allowed.
+ WriteBuffer();
+ WriteBuffer(const WriteBuffer&);
+ WriteBuffer& operator = (const WriteBuffer&);
+};
-template<class Data>
-usint
-GenericBitBuffer<Data>::reportSize()
-{
- usint bytes = sizeof(*this);
- if(this->free_buffer) { bytes += this->size * sizeof(Data); }
- return bytes;
-}
//--------------------------------------------------------------------------
file.read((char*)&(this->number_of_blocks), sizeof(this->number_of_blocks));
file.read((char*)&(this->block_size), sizeof(this->block_size));
- this->array = new usint[this->block_size * this->number_of_blocks];
- file.read((char*)(this->array), this->block_size * this->number_of_blocks * sizeof(usint));
- this->buffer = new FastBitBuffer(this->array, this->block_size);
+ usint* array_buffer = new usint[this->block_size * this->number_of_blocks];
+ file.read((char*)(array_buffer), this->block_size * this->number_of_blocks * sizeof(usint));
+ this->array = array_buffer;
this->integer_bits = length(this->size);
- this->samples = new FastBitBuffer(file, 2 * (this->number_of_blocks + 1), this->integer_bits);
+ this->samples = new ReadBuffer(file, 2 * (this->number_of_blocks + 1), this->integer_bits);
this->indexForRank();
this->indexForSelect();
std::fread(&(this->number_of_blocks), sizeof(this->number_of_blocks), 1, file);
std::fread(&(this->block_size), sizeof(this->block_size), 1, file);
- this->array = new usint[this->block_size * this->number_of_blocks];
- std::fread(this->array, sizeof(usint), this->block_size * this->number_of_blocks, file);
- this->buffer = new FastBitBuffer(this->array, this->block_size);
+ usint* array_buffer = new usint[this->block_size * this->number_of_blocks];
+ std::fread(array_buffer, sizeof(usint), this->block_size * this->number_of_blocks, file);
+ this->array = array_buffer;
this->integer_bits = length(this->size);
- this->samples = new FastBitBuffer(file, 2 * (this->number_of_blocks + 1), this->integer_bits);
+ this->samples = new ReadBuffer(file, 2 * (this->number_of_blocks + 1), this->integer_bits);
this->indexForRank();
this->indexForSelect();
}
+
BitVector::BitVector(VectorEncoder& encoder, usint universe_size) :
size(universe_size), items(encoder.items),
block_size(encoder.block_size),
number_of_blocks(encoder.blocks),
rank_index(0), select_index(0)
{
- this->array = new usint[this->block_size * this->number_of_blocks];
- this->buffer = new FastBitBuffer(this->array, this->block_size);
+ usint* array_buffer = new usint[this->block_size * this->number_of_blocks];
this->integer_bits = length(this->size);
- this->samples = new FastBitBuffer(2 * (this->number_of_blocks + 1), this->integer_bits);
+ WriteBuffer sample_buffer(2 * (this->number_of_blocks + 1), this->integer_bits);
// Copy & linearize the array.
usint pos = 0;
for(std::list<usint*>::iterator iter = encoder.array_blocks.begin(); iter != encoder.array_blocks.end(); iter++)
{
- memcpy(this->array + pos, *iter, encoder.superblock_bytes);
+ memcpy(array_buffer + pos, *iter, encoder.superblock_bytes);
pos += encoder.block_size * encoder.blocks_in_superblock;
}
- memcpy(this->array + pos, encoder.array, encoder.current_blocks * encoder.block_size * sizeof(usint));
+ memcpy(array_buffer + pos, encoder.array, encoder.current_blocks * encoder.block_size * sizeof(usint));
+ this->array = array_buffer;
// Compress the samples.
for(std::list<usint*>::iterator iter = encoder.sample_blocks.begin(); iter != encoder.sample_blocks.end(); iter++)
usint* buf = *iter;
for(usint i = 0; i < 2 * encoder.samples_in_superblock; i++)
{
- this->samples->writeItem(buf[i]);
+ sample_buffer.writeItem(buf[i]);
}
}
for(usint i = 0; i < 2 * encoder.current_samples; i++)
{
- this->samples->writeItem(encoder.samples[i]);
+ sample_buffer.writeItem(encoder.samples[i]);
}
- this->samples->writeItem(this->items);
- this->samples->writeItem(this->size);
+ sample_buffer.writeItem(this->items);
+ sample_buffer.writeItem(this->size);
+
+ this->samples = sample_buffer.getReadBuffer();
this->indexForRank();
this->indexForSelect();
BitVector::~BitVector()
{
delete[] this->array;
- delete this->buffer;
- delete this->samples;
- delete this->rank_index;
- delete this->select_index;
+ delete this->samples;
+ delete this->rank_index;
+ delete this->select_index;
}
+//--------------------------------------------------------------------------
+
void
-BitVector::writeTo(std::ofstream& file)
+BitVector::writeTo(std::ofstream& file) const
{
file.write((char*)&(this->size), sizeof(this->size));
file.write((char*)&(this->items), sizeof(this->items));
file.write((char*)&(this->number_of_blocks), sizeof(this->number_of_blocks));
file.write((char*)&(this->block_size), sizeof(this->block_size));
file.write((char*)(this->array), this->block_size * this->number_of_blocks * sizeof(usint));
- this->samples->writeBuffer(file, false);
+ this->samples->writeBuffer(file);
}
-
void
-BitVector::writeTo(FILE* file)
+BitVector::writeTo(FILE* file) const
{
std::fwrite(&(this->size), sizeof(this->size), 1, file);
std::fwrite(&(this->items), sizeof(this->items), 1, file);
std::fwrite(&(this->number_of_blocks), sizeof(this->number_of_blocks), 1, file);
std::fwrite(&(this->block_size), sizeof(this->block_size), 1, file);
std::fwrite(this->array, sizeof(usint), this->block_size * this->number_of_blocks, file);
- this->samples->writeBuffer(file, false);
+ this->samples->writeBuffer(file);
}
-
-//--------------------------------------------------------------------------
-
usint
-BitVector::reportSize()
+BitVector::reportSize() const
{
- usint bytes = this->buffer->reportSize();
- bytes += this->block_size * this->number_of_blocks * sizeof(usint);
+ // We assume the reportSize() of derived classes includes any class variables of BitVector.
+ usint bytes = this->block_size * this->number_of_blocks * sizeof(usint);
bytes += this->samples->reportSize();
bytes += this->rank_index->reportSize();
bytes += this->select_index->reportSize();
return bytes;
}
-//--------------------------------------------------------------------------
-
-usint
-BitVector::sampleForIndex(usint index)
-{
- usint low = this->select_index->readItem(index / this->select_rate);
- usint high = this->select_index->readItem(index / this->select_rate + 1);
-
- this->samples->goToItem(2 * low + 2);
- for(; low < high; low++)
- {
- if(this->samples->readItem() > index) { return low; }
- this->samples->skipItem();
- }
-
- return low;
-}
-
usint
-BitVector::sampleForValue(usint value)
+BitVector::getCompressedSize() const
{
- usint low = this->rank_index->readItem(value / this->rank_rate);
- usint high = this->rank_index->readItem(value / this->rank_rate + 1);
-
- this->samples->goToItem(2 * low + 3);
- for(; low < high; low++)
- {
- if(this->samples->readItem() > value) { return low; }
- this->samples->skipItem();
- }
-
- return low;
+ return this->block_size * this->number_of_blocks * sizeof(usint);
}
//--------------------------------------------------------------------------
usint value_samples = (this->number_of_blocks + BitVector::INDEX_RATE - 1) / BitVector::INDEX_RATE;
this->rank_rate = (this->size + value_samples - 1) / value_samples;
value_samples = (this->size + this->rank_rate - 1) / this->rank_rate + 1;
- this->rank_index = new FastBitBuffer(value_samples, length(this->number_of_blocks - 1));
+ WriteBuffer index_buffer(value_samples, length(this->number_of_blocks - 1));
usint current = 0, pointer = 0;
this->samples->goToItem(2);
usint limit = this->samples->readItem();
while(current < limit)
{
- this->rank_index->writeItem(pointer);
+ index_buffer.writeItem(pointer);
current += this->rank_rate;
}
pointer++;
}
- this->rank_index->writeItem(this->number_of_blocks - 1);
+ index_buffer.writeItem(this->number_of_blocks - 1);
+
+ this->rank_index = index_buffer.getReadBuffer();
}
void
usint index_samples = (this->number_of_blocks + BitVector::INDEX_RATE - 1) / BitVector::INDEX_RATE;
this->select_rate = (this->items + index_samples - 1) / index_samples;
index_samples = (this->items + this->select_rate - 1) / this->select_rate + 1;
- this->select_index = new FastBitBuffer(index_samples, length(this->number_of_blocks - 1));
+ WriteBuffer index_buffer(index_samples, length(this->number_of_blocks - 1));
usint current = 0, pointer = 0;
this->samples->goToItem(2);
this->samples->skipItem();
while(current < limit)
{
- this->select_index->writeItem(pointer);
+ index_buffer.writeItem(pointer);
current += this->select_rate;
}
pointer++;
}
- this->select_index->writeItem(this->number_of_blocks - 1);
+ index_buffer.writeItem(this->number_of_blocks - 1);
+
+ this->select_index = index_buffer.getReadBuffer();
+}
+
+//--------------------------------------------------------------------------
+
+BitVector::Iterator::Iterator(const BitVector& par) :
+ parent(par),
+ buffer(par.array, par.block_size),
+ samples(*(par.samples))
+{
+}
+
+BitVector::Iterator::~Iterator()
+{
+}
+
+usint
+BitVector::Iterator::sampleForIndex(usint index)
+{
+ usint low = this->parent.select_index->readItemConst(index / this->parent.select_rate);
+ usint high = this->parent.number_of_blocks - 1;
+
+ this->samples.goToItem(2 * low + 2);
+ for(; low < high; low++)
+ {
+ if(this->samples.readItem() > index) { return low; }
+ this->samples.skipItem();
+ }
+
+ return low;
+}
+
+usint
+BitVector::Iterator::sampleForValue(usint value)
+{
+ usint low = this->parent.rank_index->readItemConst(value / this->parent.rank_rate);
+ usint high = this->parent.number_of_blocks - 1;
+
+ this->samples.goToItem(2 * low + 3);
+ for(; low < high; low++)
+ {
+ if(this->samples.readItem() > value) { return low; }
+ this->samples.skipItem();
+ }
+
+ return low;
}
//--------------------------------------------------------------------------
this->samples_in_superblock = this->superblock_bytes / (2 * sizeof(usint));
this->current_samples = 0;
- this->buffer = new FastBitBuffer(this->array, this->block_size);
+ this->buffer = new WriteBuffer(this->array, this->block_size);
}
VectorEncoder::~VectorEncoder()
This should be implemented in any inherited class.
void setBit(usint value); // Values must be in increasing order.
+ void setRun(usint start, usint len);
*/
void addNewBlock();
usint size, items, blocks;
usint block_size, superblock_bytes;
- FastBitBuffer* buffer;
+ WriteBuffer* buffer;
std::list<usint*> array_blocks;
usint* array;
std::list<usint*> sample_blocks;
usint* samples;
usint samples_in_superblock, current_samples;
+
+ protected:
+
+ // These are not allowed.
+ VectorEncoder();
+ VectorEncoder(const VectorEncoder&);
+ VectorEncoder& operator = (const VectorEncoder&);
};
const static usint INDEX_RATE = 5;
BitVector(std::ifstream& file);
- BitVector(std::FILE* file);
+ BitVector(std::FILE * file);
BitVector(VectorEncoder& encoder, usint universe_size);
~BitVector();
- void writeTo(std::ofstream& file);
- void writeTo(std::FILE* file);
+//--------------------------------------------------------------------------
+
+ void writeTo(std::ofstream& file) const;
+ void writeTo(std::FILE *file) const;
+
+ inline usint getSize() const { return this->size; }
+ inline usint getNumberOfItems() const { return this->items; }
+ inline usint getBlockSize() const { return this->block_size; }
+
+ // This returns only the sizes of the dynamically allocated structures.
+ usint reportSize() const;
+
+ usint getCompressedSize() const;
//--------------------------------------------------------------------------
+ class Iterator
+ {
+ public:
+ Iterator(const BitVector& par);
+ ~Iterator();
+
+ inline bool hasNext() const
+ {
+ return (this->sample.first + this->cur < this->parent.items - 1);
+ }
+
+ protected:
+ const BitVector& parent;
+
+ usint block;
+ pair_type sample;
+ usint cur, val, run; // cur == 0 is the sample
+ usint block_items;
+
+ ReadBuffer buffer, samples;
+
+ /*
+ These functions return the sample corresponding to the block the given
+ index/value might be found in. Parameters are assumed to be valid.
+ */
+ usint sampleForIndex(usint index);
+ usint sampleForValue(usint value);
+
+ inline void getSample(usint sample_number)
+ {
+ this->block = sample_number;
+ this->samples.goToItem(2 * sample_number);
+ this->sample.first = this->samples.readItem();
+ this->sample.second = this->samples.readItem();
+ this->cur = 0;
+ this->val = this->sample.second;
+ this->block_items = this->samples.readItem() - this->sample.first - 1;
+ this->buffer.moveBuffer(this->parent.array + (this->block * this->parent.block_size));
+ }
+
+ // These are not allowed.
+ Iterator();
+ Iterator(const Iterator&);
+ Iterator& operator = (const Iterator&);
+ };
+
/*
- These should be implemented in any actual bit vector.
+ These should be implemented in any actual iterator.
- // rank is not iterator-safe
+ // rank invalidates the iterator
// regular: \sum_{i = 0}^{value} V[i]
// at_least: \sum_{i = 0}^{value - 1} V[i] + 1
usint rank(usint value, bool at_least = false);
pair_type selectRun(usint index, usint max_length);
pair_type selectNextRun(usint max_length);
- // isSet is not iterator-safe
+ // isSet invalidates the iterator
bool isSet(usint value); // V[value]
*/
- inline bool hasNext()
- {
- return (this->sample.first + cur < this->items - 1);
- }
-
//--------------------------------------------------------------------------
- inline usint getSize() { return this->size; }
- inline usint getNumberOfItems() { return this->items; }
- inline usint getBlockSize() { return this->block_size; }
-
- // This returns only the sizes of dynamically allocated structures.
- usint reportSize();
-
protected:
usint size, items;
- FastBitBuffer* buffer;
- usint* array;
- usint block_size;
- usint number_of_blocks;
+ const usint* array;
+ usint block_size;
+ usint number_of_blocks;
/*
Each sample is (rank(i) - 1, i) where V[i] = 1.
Number of samples is number_of_blocks + 1.
*/
- FastBitBuffer* samples;
- usint integer_bits;
-
- FastBitBuffer* rank_index;
- usint rank_rate;
-
- FastBitBuffer* select_index;
- usint select_rate;
+ ReadBuffer* samples;
+ usint integer_bits;
- // Iterator data. These are enough for a gap-encoded vector.
- usint block;
- pair_type sample;
- usint cur, val; // cur == 0 is the sample
- usint block_items;
+ ReadBuffer* rank_index;
+ usint rank_rate;
- /*
- These functions return the sample corresponding to the block the given
- index/value might be found in. Parameters are assumed to be valid.
- */
- usint sampleForIndex(usint index);
- usint sampleForValue(usint value);
-
- inline void getSample(usint sample_number)
- {
- this->block = sample_number;
- this->samples->goToItem(2 * sample_number);
- this->sample.first = this->samples->readItem();
- this->sample.second = this->samples->readItem();
- this->cur = 0;
- this->val = this->sample.second;
- this->block_items = this->samples->readItem() - this->sample.first - 1;
- this->buffer->moveBuffer(this->array + (this->block * this->block_size));
- }
+ ReadBuffer* select_index;
+ usint select_rate;
/*
These functions build a higher level index for faster rank/select queries.
*/
void indexForRank();
void indexForSelect();
+
+ // These are not allowed.
+ BitVector();
+ BitVector(const BitVector&);
+ BitVector& operator = (const BitVector&);
};
BitVector(file)
{
}
+
DeltaVector::DeltaVector(std::FILE * file) :
BitVector(file)
{
}
-DeltaVector::DeltaVector(DeltaEncoder& encoder, usint universe_size) :
+DeltaVector::DeltaVector(Encoder& encoder, usint universe_size) :
BitVector(encoder, universe_size)
{
}
//--------------------------------------------------------------------------
usint
-DeltaVector::rank(usint value, bool at_least)
+DeltaVector::reportSize() const
+{
+ usint bytes = sizeof(*this);
+ bytes += BitVector::reportSize();
+ return bytes;
+}
+
+//--------------------------------------------------------------------------
+
+DeltaVector::Iterator::Iterator(const DeltaVector& par) :
+ BitVector::Iterator(par)
+{
+}
+
+DeltaVector::Iterator::~Iterator()
{
- if(value >= this->size) { return this->items; }
+}
+
+usint
+DeltaVector::Iterator::rank(usint value, bool at_least)
+{
+ const DeltaVector& par = (const DeltaVector&)(this->parent);
+
+ if(value >= par.size) { return par.items; }
this->getSample(this->sampleForValue(value));
while(this->cur < this->block_items && this->val < value)
{
- this->val += this->buffer->readDeltaCode();
+ this->val += this->buffer.readDeltaCode();
this->cur++;
}
}
usint
-DeltaVector::select(usint index)
+DeltaVector::Iterator::select(usint index)
{
- if(index >= this->items) { return this->size; }
+ const DeltaVector& par = (const DeltaVector&)(this->parent);
+
+ if(index >= par.items) { return par.size; }
this->getSample(this->sampleForIndex(index));
usint lim = index - this->sample.first;
for(; this->cur < lim; this->cur++)
{
- this->val += this->buffer->readDeltaCode();
+ this->val += this->buffer.readDeltaCode();
}
return this->val;
}
usint
-DeltaVector::selectNext()
+DeltaVector::Iterator::selectNext()
{
if(this->cur >= this->block_items)
{
}
this->cur++;
- this->val += this->buffer->readDeltaCode();
+ this->val += this->buffer.readDeltaCode();
return this->val;
}
pair_type
-DeltaVector::valueAfter(usint value)
+DeltaVector::Iterator::valueAfter(usint value)
{
- if(value >= this->size) { return pair_type(this->size, this->items); }
+ const DeltaVector& par = (const DeltaVector&)(this->parent);
+
+ if(value >= par.size) { return pair_type(par.size, par.items); }
this->getSample(this->sampleForValue(value));
while(this->cur < this->block_items && this->val < value)
{
- this->val += this->buffer->readDeltaCode();
+ this->val += this->buffer.readDeltaCode();
this->cur++;
}
if(this->val < value)
}
pair_type
-DeltaVector::nextValue()
+DeltaVector::Iterator::nextValue()
{
if(this->cur >= this->block_items)
{
}
this->cur++;
- this->val += this->buffer->readDeltaCode();
+ this->val += this->buffer.readDeltaCode();
return pair_type(this->val, this->sample.first + this->cur);
}
pair_type
-DeltaVector::selectRun(usint index, usint max_length)
+DeltaVector::Iterator::selectRun(usint index, usint max_length)
{
return pair_type(this->select(index), 0);
}
pair_type
-DeltaVector::selectNextRun(usint max_length)
+DeltaVector::Iterator::selectNextRun(usint max_length)
{
return pair_type(this->selectNext(), 0);
}
bool
-DeltaVector::isSet(usint value)
+DeltaVector::Iterator::isSet(usint value)
{
- if(value >= this->size) { return false; }
+ const DeltaVector& par = (const DeltaVector&)(this->parent);
+
+ if(value >= par.size) { return false; }
this->getSample(this->sampleForValue(value));
while(this->cur < this->block_items && this->val < value)
{
- this->val += this->buffer->readDeltaCode();
+ this->val += this->buffer.readDeltaCode();
this->cur++;
}
//--------------------------------------------------------------------------
-usint
-DeltaVector::reportSize()
-{
- usint bytes = sizeof(*this);
- bytes += BitVector::reportSize();
- return bytes;
-}
-
-//--------------------------------------------------------------------------
-
DeltaEncoder::DeltaEncoder(usint block_bytes, usint superblock_size) :
VectorEncoder(block_bytes, superblock_size)
{
this->addNewBlock();
}
+void
+DeltaEncoder::setRun(usint start, usint len)
+{
+ for(usint i = start; i < start + len; i++)
+ {
+ this->setBit(i);
+ }
+}
+
+void
+DeltaEncoder::addBit(usint value)
+{
+ this->setBit(value);
+}
+
+void
+DeltaEncoder::addRun(usint start, usint len)
+{
+ this->setRun(start, len);
+}
+
+void
+DeltaEncoder::flush()
+{
+}
} // namespace CSA
~DeltaEncoder();
void setBit(usint value);
+ void setRun(usint start, usint len);
+
+ // These versions are just for compatibility with RLEVector.
+ void addBit(usint value);
+ void addRun(usint start, usint len);
+ void flush(); // Does nothing.
+
+ protected:
+
+ // These are not allowed.
+ DeltaEncoder();
+ DeltaEncoder(const DeltaEncoder&);
+ DeltaEncoder& operator = (const DeltaEncoder&);
};
class DeltaVector : public BitVector
{
public:
+ typedef DeltaEncoder Encoder;
+
DeltaVector(std::ifstream& file);
DeltaVector(std::FILE * file);
- DeltaVector(DeltaEncoder& encoder, usint universe_size);
+ DeltaVector(Encoder& encoder, usint universe_size);
~DeltaVector();
//--------------------------------------------------------------------------
- usint rank(usint value, bool at_least = false);
+ usint reportSize() const;
+
+//--------------------------------------------------------------------------
+
+ class Iterator : public BitVector::Iterator
+ {
+ public:
+ Iterator(const DeltaVector& par);
+ ~Iterator();
+
+ usint rank(usint value, bool at_least = false);
+
+ usint select(usint index);
+ usint selectNext();
- usint select(usint index);
- usint selectNext();
+ pair_type valueAfter(usint value);
+ pair_type nextValue();
- pair_type valueAfter(usint value);
- pair_type nextValue();
+ pair_type selectRun(usint index, usint max_length);
+ pair_type selectNextRun(usint max_length);
- pair_type selectRun(usint index, usint max_length);
- pair_type selectNextRun(usint max_length);
+ bool isSet(usint value);
- bool isSet(usint value);
+ // Counts the number of 1-bit runs.
+ usint countRuns();
+
+ protected:
+
+ // These are not allowed.
+ Iterator();
+ Iterator(const Iterator&);
+ Iterator& operator = (const Iterator&);
+ };
//--------------------------------------------------------------------------
- usint reportSize();
+ protected:
+
+ // These are not allowed.
+ DeltaVector();
+ DeltaVector(const DeltaVector&);
+ DeltaVector& operator = (const DeltaVector&);
};
--- /dev/null
+#include <cstdlib>
+
+#include "nibblevector.h"
+#include "../misc/utils.h"
+
+
+namespace CSA
+{
+
+
+NibbleVector::NibbleVector(std::ifstream& file) :
+ BitVector(file)
+{
+}
+
+NibbleVector::NibbleVector(Encoder& encoder, usint universe_size) :
+ BitVector(encoder, universe_size)
+{
+}
+
+NibbleVector::~NibbleVector()
+{
+}
+
+//--------------------------------------------------------------------------
+
+usint
+NibbleVector::reportSize() const
+{
+ usint bytes = sizeof(*this);
+ bytes += BitVector::reportSize();
+ return bytes;
+}
+
+//--------------------------------------------------------------------------
+
+NibbleVector::Iterator::Iterator(const NibbleVector& par) :
+ BitVector::Iterator(par),
+ use_rle(false)
+{
+}
+
+NibbleVector::Iterator::~Iterator()
+{
+}
+
+// FIXME gap encoding for all operations
+
+usint
+NibbleVector::Iterator::rank(usint value, bool at_least)
+{
+ const NibbleVector& par = (const NibbleVector&)(this->parent);
+
+ if(value >= par.size) { return par.items; }
+
+ this->valueLoop(value);
+
+ usint idx = this->sample.first + this->cur + 1;
+ if(!at_least && this->val > value)
+ {
+ idx--;
+ }
+ if(at_least && this->val < value)
+ {
+ this->getSample(this->block + 1);
+ idx = this->sample.first + this->cur + 1;
+ }
+ return idx;
+}
+
+usint
+NibbleVector::Iterator::select(usint index)
+{
+ const NibbleVector& par = (const NibbleVector&)(this->parent);
+
+ if(index >= par.items) { return par.size; }
+ this->getSample(this->sampleForIndex(index));
+
+ usint lim = index - this->sample.first;
+ while(this->cur < lim)
+ {
+ this->val += this->buffer.readNibbleCode();
+ usint temp = this->buffer.readNibbleCode();
+ this->val += temp - 1;
+ this->cur += temp;
+ }
+ if(this->cur > lim)
+ {
+ this->run = this->cur - lim;
+ this->cur -= this->run;
+ this->val -= this->run;
+ }
+
+ return this->val;
+}
+
+usint
+NibbleVector::Iterator::selectNext()
+{
+ if(this->cur >= this->block_items)
+ {
+ this->getSample(this->block + 1);
+ return this->val;
+ }
+
+ this->cur++;
+ if(this->run > 0)
+ {
+ this->val++;
+ this->run--;
+ }
+ else
+ {
+ this->val += this->buffer.readNibbleCode();
+ this->run = this->buffer.readNibbleCode() - 1;
+ }
+
+ return this->val;
+}
+
+pair_type
+NibbleVector::Iterator::valueAfter(usint value)
+{
+ const NibbleVector& par = (const NibbleVector&)(this->parent);
+
+ if(value >= par.size) { return pair_type(par.size, par.items); }
+
+ this->valueLoop(value);
+
+ if(this->val < value)
+ {
+ this->getSample(this->block + 1);
+ }
+
+ return pair_type(this->val, this->sample.first + this->cur);
+}
+
+pair_type
+NibbleVector::Iterator::nextValue()
+{
+ if(this->cur >= this->block_items)
+ {
+ this->getSample(this->block + 1);
+ return pair_type(this->val, this->sample.first);
+ }
+
+ this->cur++;
+ if(this->run > 0)
+ {
+ this->val++;
+ this->run--;
+ }
+ else
+ {
+ this->val += this->buffer.readNibbleCode();
+ this->run = this->buffer.readNibbleCode() - 1;
+ }
+
+ return pair_type(this->val, this->sample.first + this->cur);
+}
+
+pair_type
+NibbleVector::Iterator::selectRun(usint index, usint max_length)
+{
+ usint value = this->select(index);
+
+ usint len = std::min(max_length, this->run);
+ this->run -= len; this->cur += len; this->val += len;
+
+ return pair_type(value, len);
+}
+
+pair_type
+NibbleVector::Iterator::selectNextRun(usint max_length)
+{
+ usint value = this->selectNext();
+
+ usint len = std::min(max_length, this->run);
+ this->run -= len; this->cur += len; this->val += len;
+
+ return pair_type(value, len);
+}
+
+bool
+NibbleVector::Iterator::isSet(usint value)
+{
+ const NibbleVector& par = (const NibbleVector&)(this->parent);
+
+ if(value >= par.size) { return false; }
+
+ this->valueLoop(value);
+
+ return (this->val == value);
+}
+
+usint
+NibbleVector::Iterator::countRuns()
+{
+ const NibbleVector& par = (const NibbleVector&)(this->parent);
+
+ if(par.items == 0) { return 0; }
+
+ usint runs = 1;
+ pair_type res = this->selectRun(0, par.items);
+ usint last = res.first + res.second;
+
+ while(last < par.size)
+ {
+ res = this->selectNextRun(par.items);
+ if(res.first < par.size && res.first > last + 1) { runs++; }
+ last = res.first + res.second;
+ }
+
+ return runs;
+}
+
+//--------------------------------------------------------------------------
+
+NibbleEncoder::NibbleEncoder(usint block_bytes, usint superblock_size) :
+ VectorEncoder(block_bytes, superblock_size)
+{
+}
+
+NibbleEncoder::~NibbleEncoder()
+{
+}
+
+void
+NibbleEncoder::setBit(usint value)
+{
+ this->setRun(value, 1);
+}
+
+// FIXME for gap encoding
+void
+NibbleEncoder::setRun(usint start, usint len)
+{
+ if(this->items == 0)
+ {
+ this->setFirstBit(start);
+ if(len > 1)
+ {
+ this->nibbleEncode(1, len - 1);
+ }
+ return;
+ }
+ if(start < this->size || len == 0) { return; }
+
+ // Write as much into the buffer as possible.
+ usint diff = start + 1 - this->size;
+ usint free_bits = this->buffer->bitsLeft();
+ usint code_bits = this->buffer->nibbleCodeLength(diff);
+ if(free_bits >= code_bits + 4) // At least a part of the run fits into the block.
+ {
+ free_bits -= code_bits;
+ usint run_bits = this->buffer->nibbleCodeLength(len);
+ if(run_bits <= free_bits)
+ {
+ this->nibbleEncode(diff, len);
+ return;
+ }
+
+ // Encode as much as possible and let the rest spill.
+ usint llen = (usint)1 << (3 * (free_bits / 4));
+ this->nibbleEncode(diff, llen);
+ len -= llen;
+
+ // A new sample will be added.
+ this->size++;
+ this->items++;
+ }
+ else
+ {
+ this->size = start + 1;
+ this->items++;
+ }
+
+ // Didn't fit into the block. A new sample & block required.
+ this->addNewBlock();
+ if(len > 1)
+ {
+ this->nibbleEncode(1, len - 1);
+ }
+}
+
+
+} // namespace CSA
--- /dev/null
+#ifndef NIBBLEVECTOR_H
+#define NIBBLEVECTOR_H
+
+#include <fstream>
+
+#include "bitvector.h"
+
+
+namespace CSA
+{
+
+
+/*
+ This class is used to construct a NibbleVector.
+*/
+
+class NibbleEncoder : public VectorEncoder
+{
+ public:
+ NibbleEncoder(usint block_bytes, usint superblock_size = VectorEncoder::SUPERBLOCK_SIZE);
+ ~NibbleEncoder();
+
+ void setBit(usint value);
+ void setRun(usint start, usint len);
+
+ // FIXME for gap encoding
+ inline void nibbleEncode(usint diff, usint len)
+ {
+ this->size += diff + len - 1;
+ this->items += len;
+ this->buffer->writeNibbleCode(diff);
+ this->buffer->writeNibbleCode(len);
+ }
+
+ protected:
+
+ // These are not allowed.
+ NibbleEncoder();
+ NibbleEncoder(const NibbleEncoder&);
+ NibbleEncoder& operator = (const NibbleEncoder&);
+};
+
+
+/*
+ This bit vector uses nibble coding. Each block is either run-length encoded or
+ gap encoded, depending on the first nibble.
+
+ // FIXME reverting to gap encoding not implemented yet
+*/
+
+class NibbleVector : public BitVector
+{
+ public:
+ typedef NibbleEncoder Encoder;
+
+ NibbleVector(std::ifstream& file);
+ NibbleVector(Encoder& encoder, usint universe_size);
+ ~NibbleVector();
+
+//--------------------------------------------------------------------------
+
+ usint reportSize() const;
+
+//--------------------------------------------------------------------------
+
+ class Iterator : public BitVector::Iterator
+ {
+ public:
+ Iterator(const NibbleVector& par);
+ ~Iterator();
+
+ usint rank(usint value, bool at_least = false);
+
+ usint select(usint index);
+ usint selectNext();
+
+ pair_type valueAfter(usint value);
+ pair_type nextValue();
+
+ pair_type selectRun(usint index, usint max_length);
+ pair_type selectNextRun(usint max_length);
+
+ bool isSet(usint value);
+
+ // Counts the number of 1-bit runs.
+ usint countRuns();
+
+ static const usint GAP_ENCODING = 0;
+ static const usint RUN_LENGTH_ENCODING = 1;
+
+ protected:
+
+ // FIXME for gap encoding
+ inline void valueLoop(usint value)
+ {
+ this->getSample(this->sampleForValue(value));
+
+ if(this->val >= value) { return; }
+ while(this->cur < this->block_items)
+ {
+ this->val += this->buffer.readNibbleCode();
+ this->cur++;
+ this->run = this->buffer.readNibbleCode() - 1;
+ if(this->val >= value) { break; }
+
+ this->cur += this->run;
+ this->val += this->run;
+ if(this->val >= value)
+ {
+ this->run = this->val - value;
+ this->val = value;
+ this->cur -= this->run;
+ break;
+ }
+ this->run = 0;
+ }
+ }
+
+ inline void getSample(usint sample_number)
+ {
+ BitVector::Iterator::getSample(sample_number);
+ this->run = 0;
+// this->use_rle = this->buffer.readNibble();
+ }
+
+ bool use_rle;
+
+ // These are not allowed.
+ Iterator();
+ Iterator(const Iterator&);
+ Iterator& operator = (const Iterator&);
+ };
+
+//--------------------------------------------------------------------------
+
+ protected:
+
+ // These are not allowed.
+ NibbleVector();
+ NibbleVector(const NibbleVector&);
+ NibbleVector& operator = (const NibbleVector&);
+};
+
+
+} // namespace CSA
+
+
+#endif // NIBBLEVECTOR_H
{
}
-RLEVector::RLEVector(RLEEncoder& encoder, usint universe_size) :
+RLEVector::RLEVector(Encoder& encoder, usint universe_size) :
BitVector(encoder, universe_size)
{
}
//--------------------------------------------------------------------------
usint
-RLEVector::rank(usint value, bool at_least)
+RLEVector::reportSize() const
{
- if(value >= this->size) { return this->items; }
+ usint bytes = sizeof(*this);
+ bytes += BitVector::reportSize();
+ return bytes;
+}
+
+//--------------------------------------------------------------------------
+
+RLEVector::Iterator::Iterator(const RLEVector& par) :
+ BitVector::Iterator(par)
+{
+}
+
+RLEVector::Iterator::~Iterator()
+{
+}
+
+usint
+RLEVector::Iterator::rank(usint value, bool at_least)
+{
+ const RLEVector& par = (const RLEVector&)(this->parent);
+
+ if(value >= par.size) { return par.items; }
this->valueLoop(value);
}
usint
-RLEVector::select(usint index)
+RLEVector::Iterator::select(usint index)
{
- if(index >= this->items) { return this->size; }
+ const RLEVector& par = (const RLEVector&)(this->parent);
+
+ if(index >= par.items) { return par.size; }
this->getSample(this->sampleForIndex(index));
this->run = 0;
usint lim = index - this->sample.first;
while(this->cur < lim)
{
- this->val += this->buffer->readDeltaCode();
- usint temp = this->buffer->readDeltaCode();
+ this->val += this->buffer.readDeltaCode();
+ usint temp = this->buffer.readDeltaCode();
this->val += temp - 1;
this->cur += temp;
}
}
usint
-RLEVector::selectNext()
+RLEVector::Iterator::selectNext()
{
if(this->cur >= this->block_items)
{
}
else
{
- this->val += this->buffer->readDeltaCode();
- this->run = this->buffer->readDeltaCode() - 1;
+ this->val += this->buffer.readDeltaCode();
+ this->run = this->buffer.readDeltaCode() - 1;
}
return this->val;
}
pair_type
-RLEVector::valueAfter(usint value)
+RLEVector::Iterator::valueAfter(usint value)
{
- if(value >= this->size) { return pair_type(this->size, this->items); }
+ const RLEVector& par = (const RLEVector&)(this->parent);
+
+ if(value >= par.size) { return pair_type(par.size, par.items); }
this->valueLoop(value);
}
pair_type
-RLEVector::nextValue()
+RLEVector::Iterator::nextValue()
{
if(this->cur >= this->block_items)
{
}
else
{
- this->val += this->buffer->readDeltaCode();
- this->run = this->buffer->readDeltaCode() - 1;
+ this->val += this->buffer.readDeltaCode();
+ this->run = this->buffer.readDeltaCode() - 1;
}
return pair_type(this->val, this->sample.first + this->cur);
}
pair_type
-RLEVector::selectRun(usint index, usint max_length)
+RLEVector::Iterator::selectRun(usint index, usint max_length)
{
usint value = this->select(index);
}
pair_type
-RLEVector::selectNextRun(usint max_length)
+RLEVector::Iterator::selectNextRun(usint max_length)
{
usint value = this->selectNext();
}
bool
-RLEVector::isSet(usint value)
+RLEVector::Iterator::isSet(usint value)
{
- if(value >= this->size) { return false; }
+ const RLEVector& par = (const RLEVector&)(this->parent);
+
+ if(value >= par.size) { return false; }
this->valueLoop(value);
return (this->val == value);
}
-//--------------------------------------------------------------------------
-
usint
-RLEVector::reportSize()
+RLEVector::Iterator::countRuns()
{
- usint bytes = sizeof(*this);
- bytes += BitVector::reportSize();
- return bytes;
+ const RLEVector& par = (const RLEVector&)(this->parent);
+
+ if(par.items == 0) { return 0; }
+
+ usint runs = 1;
+ pair_type res = this->selectRun(0, par.items);
+ usint last = res.first + res.second;
+
+ while(last < par.size)
+ {
+ res = this->selectNextRun(par.items);
+ if(res.first < par.size && res.first > last + 1) { runs++; }
+ last = res.first + res.second;
+ }
+
+ return runs;
}
//--------------------------------------------------------------------------
RLEEncoder::RLEEncoder(usint block_bytes, usint superblock_size) :
- VectorEncoder(block_bytes, superblock_size)
+ VectorEncoder(block_bytes, superblock_size),
+ run(EMPTY_PAIR)
{
}
{
}
+void
+RLEEncoder::setBit(usint value)
+{
+ this->setRun(value, 1);
+}
+
void
RLEEncoder::setRun(usint start, usint len)
{
}
}
+void
+RLEEncoder::addBit(usint value)
+{
+ this->addRun(value, 1);
+}
+
+void
+RLEEncoder::addRun(usint start, usint len)
+{
+ if(this->run.second == 0)
+ {
+ this->run = pair_type(start, len);
+ }
+ else if(start == this->run.first + this->run.second)
+ {
+ this->run.second += len;
+ }
+ else
+ {
+ this->setRun(this->run.first, this->run.second);
+ this->run = pair_type(start, len);
+ }
+}
+
+void
+RLEEncoder::flush()
+{
+ this->setRun(this->run.first, this->run.second);
+ this->run.second = 0;
+}
+
} // namespace CSA
RLEEncoder(usint block_bytes, usint superblock_size = VectorEncoder::SUPERBLOCK_SIZE);
~RLEEncoder();
-// void setBit(usint value);
+ void setBit(usint value);
void setRun(usint start, usint len);
+ // These versions try to combine the runs if possible.
+ void addBit(usint value);
+ void addRun(usint start, usint len);
+ void flush(); // Call this when finished.
+
inline void RLEncode(usint diff, usint len)
{
this->size += diff + len - 1;
this->buffer->writeDeltaCode(diff);
this->buffer->writeDeltaCode(len);
}
+
+ protected:
+ pair_type run;
+
+ // These are not allowed.
+ RLEEncoder();
+ RLEEncoder(const RLEEncoder&);
+ RLEEncoder& operator = (const RLEEncoder&);
};
class RLEVector : public BitVector
{
public:
+ typedef RLEEncoder Encoder;
+
RLEVector(std::ifstream& file);
- RLEVector(RLEEncoder& encoder, usint universe_size);
+ RLEVector(Encoder& encoder, usint universe_size);
~RLEVector();
//--------------------------------------------------------------------------
- usint rank(usint value, bool at_least = false);
+ usint reportSize() const;
- usint select(usint index);
- usint selectNext();
+//--------------------------------------------------------------------------
- pair_type valueAfter(usint value);
- pair_type nextValue();
+ class Iterator : public BitVector::Iterator
+ {
+ public:
+ Iterator(const RLEVector& par);
+ ~Iterator();
- pair_type selectRun(usint index, usint max_length);
- pair_type selectNextRun(usint max_length);
+ usint rank(usint value, bool at_least = false);
- bool isSet(usint value);
+ usint select(usint index);
+ usint selectNext();
-//--------------------------------------------------------------------------
+ pair_type valueAfter(usint value);
+ pair_type nextValue();
- usint reportSize();
+ pair_type selectRun(usint index, usint max_length);
+ pair_type selectNextRun(usint max_length);
- protected:
- usint run;
+ bool isSet(usint value);
- inline void valueLoop(usint value)
- {
- this->getSample(this->sampleForValue(value));
- this->run = 0;
-
- if(this->val >= value) { return; }
- while(this->cur < this->block_items)
- {
- this->val += this->buffer->readDeltaCode();
- this->cur++;
- this->run = this->buffer->readDeltaCode() - 1;
- if(this->val >= value) { break; }
-
- this->cur += this->run;
- this->val += this->run;
- if(this->val >= value)
+ // Counts the number of 1-bit runs.
+ usint countRuns();
+
+ protected:
+
+ inline void valueLoop(usint value)
{
- this->run = this->val - value;
- this->val = value;
- this->cur -= this->run;
- break;
+ this->getSample(this->sampleForValue(value));
+ this->run = 0;
+
+ if(this->val >= value) { return; }
+ while(this->cur < this->block_items)
+ {
+ this->val += this->buffer.readDeltaCode();
+ this->cur++;
+ this->run = this->buffer.readDeltaCode() - 1;
+ if(this->val >= value) { break; }
+
+ this->cur += this->run;
+ this->val += this->run;
+ if(this->val >= value)
+ {
+ this->run = this->val - value;
+ this->val = value;
+ this->cur -= this->run;
+ break;
+ }
+ this->run = 0;
+ }
}
- this->run = 0;
- }
- }
+
+ // These are not allowed.
+ Iterator();
+ Iterator(const Iterator&);
+ Iterator& operator = (const Iterator&);
+ };
+
+//--------------------------------------------------------------------------
+
+ protected:
+
+ // These are not allowed.
+ RLEVector();
+ RLEVector(const RLEVector&);
+ RLEVector& operator = (const RLEVector&);
};
{
if((first == 0 && second == 0) || positions == 0) { return 0; }
+ RLEVector::Iterator* first_iter = 0;
+ RLEVector::Iterator* second_iter = 0;
+
pair_type first_run;
bool first_finished;
if(first == 0)
}
else
{
- first_run = first->selectRun(0, size);
+ first_iter = new RLEVector::Iterator(*first);
+ first_run = first_iter->selectRun(0, size);
first_run.second++;
first_finished = false;
}
}
else
{
- second_bit = second->select(0);
+ second_iter = new RLEVector::Iterator(*second);
+ second_bit = second_iter->select(0);
}
RLEEncoder encoder(block_size);
handleRun(encoder, run, first_run, positions[i]);
if(first_run.second == 0)
{
- if(first->hasNext())
+ if(first_iter->hasNext())
{
- first_run = first->selectNextRun(size);
+ first_run = first_iter->selectNextRun(size);
first_run.first += i;
first_run.second++;
}
if(i == second_bit) // positions[i] is one
{
handleOne(encoder, run, positions[i]);
- second_bit = second->selectNext();
+ second_bit = second_iter->selectNext();
}
else // positions[i] is zero
{
while(!first_finished)
{
handleRun(encoder, run, first_run, size);
- if(first->hasNext())
+ if(first_iter->hasNext())
{
- first_run = first->selectNextRun(size);
+ first_run = first_iter->selectNextRun(size);
first_run.first += n;
first_run.second++;
}
encoder.setRun(run.first, run.second);
}
+ delete first_iter; delete second_iter;
delete first; delete second;
return new RLEVector(encoder, size);
}
{
if((first == 0 && second == 0) || positions == 0) { return 0; }
+ DeltaVector::Iterator* first_iter = 0;
+ DeltaVector::Iterator* second_iter = 0;
+
usint first_bit;
bool first_finished;
if(first == 0)
}
else
{
- first_bit = first->select(0);
+ first_iter = new DeltaVector::Iterator(*first);
+ first_bit = first_iter->select(0);
first_finished = false;
}
}
else
{
- second_bit = second->select(0);
+ second_iter = new DeltaVector::Iterator(*second);
+ second_bit = second_iter->select(0);
}
DeltaEncoder encoder(block_size);
while(!first_finished && first_bit < positions[i])
{
encoder.setBit(first_bit);
- if(first->hasNext())
+ if(first_iter->hasNext())
{
- first_bit = first->selectNext() + i;
+ first_bit = first_iter->selectNext() + i;
}
else
{
if(i == second_bit) // positions[i] is one
{
encoder.setBit(positions[i]);
- second_bit = second->selectNext();
+ second_bit = second_iter->selectNext();
}
}
while(!first_finished)
{
encoder.setBit(first_bit);
- if(!first->hasNext()) { break; }
- first_bit = first->selectNext() + n;
+ if(!first_iter->hasNext()) { break; }
+ first_bit = first_iter->selectNext() + n;
}
+ delete first_iter; delete second_iter;
delete first; delete second;
return new DeltaVector(encoder, size);
}
#ifndef VECTORS_H
#define VECTORS_H
-#include "deltavector.h"
-#include "rlevector.h"
-
namespace CSA
{
The original vectors are deleted.
*/
-RLEVector* mergeVectors(RLEVector* first, RLEVector* second, usint* positions, usint n, usint size, usint block_size);
+template<class V, class E>
+void
+handleOne(E& encoder, pair_type& run, usint position)
+{
+ if(run.second == 0)
+ {
+ run.first = position;
+ run.second = 1;
+ return;
+ }
+ if(position == run.first + run.second)
+ {
+ run.second++;
+ return;
+ }
+ encoder.setRun(run.first, run.second);
+ run.first = position;
+ run.second = 1;
+}
+
+
+template<class V, class E>
+void
+handleRun(E& encoder, pair_type& run, pair_type& next, usint limit)
+{
+ if(run.second == 0)
+ {
+ run.first = next.first;
+ run.second = std::min(limit - run.first, next.second);
+ next.first += run.second;
+ next.second -= run.second;
+ return;
+ }
+
+ if(next.first == run.first + run.second)
+ {
+ usint cont = std::min(limit - next.first, next.second);
+ run.second += cont;
+ next.first += cont;
+ next.second -= cont;
+ return;
+ }
+
+ encoder.setRun(run.first, run.second);
+ run.first = next.first;
+ run.second = std::min(limit - run.first, next.second);;
+ next.first += run.second;
+ next.second -= run.second;
+}
+
+
+template<class V, class E, class I>
+V*
+mergeVectors(V* first, V* second, usint* positions, usint n, usint size, usint block_size)
+{
+ if((first == 0 && second == 0) || positions == 0) { return 0; }
+
+ I* first_iter = 0;
+ I* second_iter = 0;
+
+ pair_type first_run;
+ bool first_finished;
+ if(first == 0)
+ {
+ first_run = pair_type(size, 0);
+ first_finished = true;
+ }
+ else
+ {
+ first_iter = new I(*first);
+ first_run = first_iter->selectRun(0, size);
+ first_run.second++;
+ first_finished = false;
+ }
+
+ usint second_bit;
+ if(second == 0)
+ {
+ second_bit = n;
+ }
+ else
+ {
+ second_iter = new I(*second);
+ second_bit = second_iter->select(0);
+ }
+
+ E encoder(block_size);
+ pair_type run = pair_type(size, 0);
+ for(usint i = 0; i < n; i++, first_run.first++)
+ {
+ while(!first_finished && first_run.first < positions[i])
+ {
+ handleRun<V, E>(encoder, run, first_run, positions[i]);
+ if(first_run.second == 0)
+ {
+ if(first_iter->hasNext())
+ {
+ first_run = first_iter->selectNextRun(size);
+ first_run.first += i;
+ first_run.second++;
+ }
+ else
+ {
+ first_finished = true;
+ }
+ }
+ }
+
+ if(i == second_bit) // positions[i] is one
+ {
+ handleOne<V, E>(encoder, run, positions[i]);
+ second_bit = second_iter->selectNext();
+ }
+ else // positions[i] is zero
+ {
+ if(run.second != 0)
+ {
+ encoder.setRun(run.first, run.second);
+ run.second = 0;
+ }
+ }
+ }
+
+ while(!first_finished)
+ {
+ handleRun<V, E>(encoder, run, first_run, size);
+ if(first_iter->hasNext())
+ {
+ first_run = first_iter->selectNextRun(size);
+ first_run.first += n;
+ first_run.second++;
+ }
+ else { break; }
+ }
+
+ if(run.second != 0)
+ {
+ encoder.setRun(run.first, run.second);
+ }
-DeltaVector* mergeVectors(DeltaVector* first, DeltaVector* second, usint* positions, usint n, usint size, usint block_size);
+ delete first_iter; delete second_iter;
+ delete first; delete second;
+ return new V(encoder, size);
+}
} // namespace CSA
--- /dev/null
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+
+
+using namespace CSA;
+
+
+/*
+ This program writes run-length encoded PLCP of the collection into a file.
+*/
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "PLCP builder" << std::endl;
+ if(argc < 2)
+ {
+ std::cout << "Usage: build_plcp base_name [block_size]" << std::endl;
+ return 1;
+ }
+
+ std::string base_name = argv[1];
+ std::string plcp_name = base_name + PLCP_EXTENSION;
+ std::cout << "PLCP: " << plcp_name << std::endl;
+ std::ofstream plcp_file(plcp_name.c_str(), std::ios_base::binary);
+ if(!plcp_file)
+ {
+ std::cerr << "Error creating PLCP file!" << std::endl;
+ return 2;
+ }
+
+ usint block_size = 32;
+ if(argc > 2) { block_size = atoi(argv[2]); }
+ std::cout << "Block size: " << block_size << std::endl;
+ std::cout << std::endl;
+ RLCSA rlcsa(base_name);
+
+ clock_t start = clock();
+ RLEVector* plcp = rlcsa.buildPLCP(block_size);
+ plcp->writeTo(plcp_file);
+ clock_t stop = clock();
+
+ double time = ((stop - start) / (double)CLOCKS_PER_SEC);
+ double megabytes = rlcsa.getSize() / (double)MEGABYTE;
+ double size = plcp->reportSize() / (double)MEGABYTE;
+ RLEVector::Iterator iter(*plcp);
+ usint runs = iter.countRuns();
+
+ std::cout << megabytes << " megabytes in " << time << " seconds (" << (megabytes / time) << " MB/s)" << std::endl;
+ std::cout << "PLCP size: " << size << " MB" << std::endl;
+ std::cout << "Runs: " << runs << std::endl;
+ std::cout << std::endl;
+
+ plcp_file.close();
+ delete plcp;
+
+ return 0;
+}
--- /dev/null
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+#include "rlcsa_builder.h"
+#include "misc/utils.h"
+
+
+using namespace CSA;
+
+
+
+int
+lineByLineRLCSA(std::string base_name, usint block_size, usint sample_rate, usint buffer_size)
+{
+ Parameters parameters;
+ parameters.set(RLCSA_BLOCK_SIZE.first, block_size);
+ parameters.set(SAMPLE_RATE.first, sample_rate);
+
+ if(sample_rate > 0)
+ {
+ parameters.set(SUPPORT_LOCATE.first, 1);
+ parameters.set(SUPPORT_DISPLAY.first, 1);
+ }
+ else
+ {
+ parameters.set(SUPPORT_LOCATE.first, 0);
+ parameters.set(SUPPORT_DISPLAY.first, 0);
+ }
+
+ std::string parameters_name = base_name + PARAMETERS_EXTENSION;
+ parameters.print();
+ parameters.write(parameters_name);
+
+ std::cout << "Input: " << base_name << std::endl;
+ std::ifstream input_file(base_name.c_str(), std::ios_base::binary);
+ if(!input_file)
+ {
+ std::cerr << "Error opening input file!" << std::endl;
+ return 2;
+ }
+ std::cout << "Buffer size: " << buffer_size << " MB" << std::endl;
+ std::cout << std::endl;
+
+ double start = readTimer();
+ RLCSABuilder builder(parameters.get(RLCSA_BLOCK_SIZE), parameters.get(SAMPLE_RATE), buffer_size * MEGABYTE);
+
+ usint lines = 0, total = 0;
+ while(input_file)
+ {
+ char buffer[16384]; // FIXME What if lines are longer? Probably fails.
+ input_file.getline(buffer, 16384);
+ usint chars = input_file.gcount();
+ lines++; total += chars;
+ if(chars >= 16383) { std::cout << lines << ": " << chars << " chars read!" << std::endl; }
+ if(chars > 1) { builder.insertSequence(buffer, chars - 1, false); }
+ }
+
+ RLCSA* rlcsa = 0;
+ if(builder.isOk())
+ {
+ rlcsa = builder.getRLCSA();
+ rlcsa->writeTo(base_name);
+ }
+ else
+ {
+ std::cerr << "Error: RLCSA construction failed!" << std::endl;
+ return 3;
+ }
+
+ double time = readTimer() - start;
+ double build_time = builder.getBuildTime();
+ double search_time = builder.getSearchTime();
+ double sort_time = builder.getSortTime();
+ double merge_time = builder.getMergeTime();
+
+ double megabytes = rlcsa->getSize() / (double)MEGABYTE;
+ usint sequences = rlcsa->getNumberOfSequences();
+ std::cout << sequences << " sequences" << std::endl;
+ std::cout << megabytes << " megabytes in " << time << " seconds (" << (megabytes / time) << " MB/s)" << std::endl;
+ std::cout << "(build " << build_time << " s, search " << search_time << "s, sort " << sort_time << " s, merge " << merge_time << " s)" << std::endl;
+ std::cout << std::endl;
+
+ delete rlcsa;
+ return 0;
+}
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "Line-by-line RLCSA builder" << std::endl;
+ if(argc < 3)
+ {
+ std::cout << "Usage: build_rlcsa base_name buffer_size [block_size [sample_rate]]" << std::endl;
+ return 1;
+ }
+ std::cout << std::endl;
+
+ int name_arg = 1, buffer_arg = 2, block_arg = 3, sample_arg = 4;
+ usint block_size = (argc > block_arg ? atoi(argv[block_arg]) : RLCSA_BLOCK_SIZE.second);
+ usint sample_rate = (argc > sample_arg ? atoi(argv[sample_arg]) : SAMPLE_RATE.second);
+ return lineByLineRLCSA(argv[name_arg], block_size, sample_rate, atoi(argv[buffer_arg]));
+}
-rlcsa.o: rlcsa.cpp rlcsa.h bits/vectors.h bits/deltavector.h \
+build_plcp.o: build_plcp.cpp rlcsa.h bits/deltavector.h bits/bitvector.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/rlevector.h \
+ bits/nibblevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
+ lcpsamples.h bits/array.h bits/vectors.h misc/utils.h misc/definitions.h \
+ misc/parameters.h
+build_rlcsa.o: build_rlcsa.cpp rlcsa_builder.h rlcsa.h bits/deltavector.h \
bits/bitvector.h bits/../misc/definitions.h bits/bitbuffer.h \
- bits/rlevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
- bits/deltavector.h misc/parameters.h misc/definitions.h misc/utils.h \
- qsufsort/qsufsort.h qsufsort/../misc/definitions.h
-rlcsa_builder.o: rlcsa_builder.cpp rlcsa_builder.h rlcsa.h bits/vectors.h \
+ bits/rlevector.h bits/nibblevector.h sasamples.h misc/definitions.h \
+ bits/bitbuffer.h lcpsamples.h bits/array.h bits/vectors.h misc/utils.h \
+ misc/definitions.h misc/parameters.h
+display_test.o: display_test.cpp rlcsa.h bits/deltavector.h \
+ bits/bitvector.h bits/../misc/definitions.h bits/bitbuffer.h \
+ bits/rlevector.h bits/nibblevector.h sasamples.h misc/definitions.h \
+ bits/bitbuffer.h lcpsamples.h bits/array.h bits/vectors.h misc/utils.h \
+ misc/definitions.h misc/parameters.h
+extract_sequence.o: extract_sequence.cpp rlcsa.h bits/deltavector.h \
+ bits/bitvector.h bits/../misc/definitions.h bits/bitbuffer.h \
+ bits/rlevector.h bits/nibblevector.h sasamples.h misc/definitions.h \
+ bits/bitbuffer.h lcpsamples.h bits/array.h bits/vectors.h misc/utils.h \
+ misc/definitions.h misc/parameters.h
+lcp_test.o: lcp_test.cpp rlcsa.h bits/deltavector.h bits/bitvector.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/rlevector.h \
+ bits/nibblevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
+ lcpsamples.h bits/array.h bits/vectors.h misc/utils.h misc/definitions.h \
+ misc/parameters.h
+lcpsamples.o: lcpsamples.cpp lcpsamples.h bits/array.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/vectors.h misc/utils.h \
+ misc/definitions.h
+locate_test.o: locate_test.cpp rlcsa.h bits/deltavector.h \
+ bits/bitvector.h bits/../misc/definitions.h bits/bitbuffer.h \
+ bits/rlevector.h bits/nibblevector.h sasamples.h misc/definitions.h \
+ bits/bitbuffer.h lcpsamples.h bits/array.h bits/vectors.h misc/utils.h \
+ misc/definitions.h misc/parameters.h
+parallel_build.o: parallel_build.cpp rlcsa_builder.h rlcsa.h \
+ bits/deltavector.h bits/bitvector.h bits/../misc/definitions.h \
+ bits/bitbuffer.h bits/rlevector.h bits/nibblevector.h sasamples.h \
+ misc/definitions.h bits/bitbuffer.h lcpsamples.h bits/array.h \
+ bits/vectors.h misc/utils.h misc/definitions.h misc/parameters.h
+read_bwt.o: read_bwt.cpp rlcsa.h bits/deltavector.h bits/bitvector.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/rlevector.h \
+ bits/nibblevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
+ lcpsamples.h bits/array.h bits/vectors.h misc/utils.h misc/definitions.h \
+ misc/parameters.h
+rlcsa.o: rlcsa.cpp rlcsa.h bits/deltavector.h bits/bitvector.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/rlevector.h \
+ bits/nibblevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
+ lcpsamples.h bits/array.h bits/vectors.h misc/utils.h misc/definitions.h \
+ misc/parameters.h qsufsort/qsufsort.h qsufsort/../misc/definitions.h
+rlcsa_builder.o: rlcsa_builder.cpp rlcsa_builder.h rlcsa.h \
bits/deltavector.h bits/bitvector.h bits/../misc/definitions.h \
- bits/bitbuffer.h bits/rlevector.h sasamples.h misc/definitions.h \
- bits/bitbuffer.h bits/deltavector.h misc/parameters.h misc/definitions.h
+ bits/bitbuffer.h bits/rlevector.h bits/nibblevector.h sasamples.h \
+ misc/definitions.h bits/bitbuffer.h lcpsamples.h bits/array.h \
+ bits/vectors.h misc/utils.h misc/definitions.h misc/parameters.h
+rlcsa_grep.o: rlcsa_grep.cpp rlcsa.h bits/deltavector.h bits/bitvector.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/rlevector.h \
+ bits/nibblevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
+ lcpsamples.h bits/array.h bits/vectors.h misc/utils.h misc/definitions.h \
+ misc/parameters.h
+rlcsa_test.o: rlcsa_test.cpp rlcsa.h bits/deltavector.h bits/bitvector.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/rlevector.h \
+ bits/nibblevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
+ lcpsamples.h bits/array.h bits/vectors.h misc/utils.h misc/definitions.h \
+ misc/parameters.h
+sample_lcp.o: sample_lcp.cpp rlcsa.h bits/deltavector.h bits/bitvector.h \
+ bits/../misc/definitions.h bits/bitbuffer.h bits/rlevector.h \
+ bits/nibblevector.h sasamples.h misc/definitions.h bits/bitbuffer.h \
+ lcpsamples.h bits/array.h bits/vectors.h misc/utils.h misc/definitions.h \
+ misc/parameters.h
sasamples.o: sasamples.cpp sasamples.h misc/definitions.h \
bits/bitbuffer.h bits/../misc/definitions.h bits/deltavector.h \
bits/bitvector.h bits/bitbuffer.h misc/utils.h misc/definitions.h
+array.o: bits/array.cpp bits/array.h bits/../misc/definitions.h \
+ bits/bitbuffer.h
bitvector.o: bits/bitvector.cpp bits/bitvector.h \
bits/../misc/definitions.h bits/bitbuffer.h
deltavector.o: bits/deltavector.cpp bits/deltavector.h bits/bitvector.h \
bits/../misc/definitions.h bits/bitbuffer.h
+nibblevector.o: bits/nibblevector.cpp bits/nibblevector.h \
+ bits/bitvector.h bits/../misc/definitions.h bits/bitbuffer.h \
+ bits/../misc/utils.h bits/../misc/definitions.h
rlevector.o: bits/rlevector.cpp bits/rlevector.h bits/bitvector.h \
bits/../misc/definitions.h bits/bitbuffer.h bits/../misc/utils.h \
bits/../misc/definitions.h
-vectors.o: bits/vectors.cpp bits/vectors.h bits/deltavector.h \
- bits/bitvector.h bits/../misc/definitions.h bits/bitbuffer.h \
- bits/rlevector.h bits/../misc/utils.h bits/../misc/definitions.h
+vectors.o: bits/vectors.cpp bits/vectors.h bits/../misc/utils.h \
+ bits/../misc/definitions.h
parameters.o: misc/parameters.cpp misc/parameters.h misc/definitions.h
utils.o: misc/utils.cpp misc/utils.h misc/definitions.h
qsufsort.o: qsufsort/qsufsort.c qsufsort/qsufsort.h \
qsufsort/../misc/definitions.h
+convert_patterns.o: utils/convert_patterns.cpp \
+ utils/../misc/definitions.h utils/../misc/utils.h \
+ utils/../misc/definitions.h
+extract_text.o: utils/extract_text.cpp utils/../misc/definitions.h
+split_text.o: utils/split_text.cpp utils/../misc/definitions.h \
+ utils/../misc/utils.h utils/../misc/definitions.h
--- /dev/null
+#include <algorithm>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+#ifdef MULTITHREAD_SUPPORT
+#include <omp.h>
+const int MAX_THREADS = 64;
+#endif
+
+
+using namespace CSA;
+
+
+int main(int argc, char** argv)
+{
+ std::cout << "RLCSA display test" << std::endl;
+ if(argc < 4)
+ {
+ std::cout << "Usage: display_test basename sequences max_length [threads [random_seed]]" << std::endl;
+ return 1;
+ }
+
+ std::cout << "Base name: " << argv[1] << std::endl;
+
+ sint sequences = std::max(atoi(argv[2]), 1);
+ std::cout << "Sequences: " << sequences << std::endl;
+
+ usint max_length = std::max(atoi(argv[3]), 1);
+ std::cout << "Prefix length: " << max_length << std::endl;
+
+ sint threads = 1;
+#ifdef MULTITHREAD_SUPPORT
+ if(argc > 4)
+ {
+ threads = std::min(MAX_THREADS, std::max(atoi(argv[4]), 1));
+ }
+#endif
+ std::cout << "Threads: " << threads << std::endl;
+
+ usint seed = 0xDEADBEEF;
+ if(argc > 5)
+ {
+ seed = atoi(argv[5]);
+ }
+ std::cout << "Random seed: " << seed << std::endl;
+ std::cout << std::endl;
+
+ RLCSA rlcsa(argv[1]);
+ if(!rlcsa.supportsDisplay())
+ {
+ std::cerr << "Error: Display is not supported!" << std::endl;
+ return 2;
+ }
+ rlcsa.printInfo();
+ rlcsa.reportSize(true);
+
+ usint total = 0;
+ usint seq_num, seq_total = rlcsa.getNumberOfSequences();
+ sint i;
+
+ double start = readTimer();
+ srand(seed);
+ uchar* buffer = new uchar[max_length];
+ #ifdef MULTITHREAD_SUPPORT
+ usint length;
+ usint thread_id;
+ uchar* buffers[threads]; buffers[0] = buffer;
+ for(i = 1; i < threads; i++)
+ {
+ buffers[i] = new uchar[max_length];
+ }
+ omp_set_num_threads(threads);
+ #pragma omp parallel private(seq_num, length, thread_id)
+ {
+ #pragma omp for schedule(dynamic, 1)
+ for(i = 0; i < sequences; i++)
+ {
+ #pragma omp critical
+ {
+ seq_num = rand() % seq_total;
+ }
+ thread_id = omp_get_thread_num();
+ length = rlcsa.displayPrefix(seq_num, max_length, buffers[thread_id]);
+ #pragma omp critical
+ {
+ total += length;
+ }
+ }
+ }
+ for(i = 1; i < threads; i++) { delete[] buffers[i]; }
+ #else
+ for(i = 0; i < sequences; i++)
+ {
+ seq_num = rand() % seq_total;
+ total += rlcsa.displayPrefix(seq_num, max_length, buffer);
+ }
+ #endif
+ delete[] buffer;
+
+ double time = readTimer() - start;
+ double megabytes = total / (double)MEGABYTE;
+ std::cout << megabytes << " megabytes in " << time << " seconds (" << (megabytes / time) << " MB/s)" << std::endl;
+ std::cout << std::endl;
+
+ return 0;
+}
--- /dev/null
+#include <algorithm>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+#ifdef MULTITHREAD_SUPPORT
+#include <omp.h>
+const int MAX_THREADS = 64;
+#endif
+
+
+using namespace CSA;
+
+
+int main(int argc, char** argv)
+{
+ std::cout << "RLCSA display test" << std::endl;
+ if(argc < 4)
+ {
+ std::cout << "Usage: extract_sequence base_name sequence_number output" << std::endl;
+ return 1;
+ }
+
+ std::cout << "Base name: " << argv[1] << std::endl;
+ usint sequence = atoi(argv[2]);
+ std::cout << "Sequence number: " << sequence << std::endl;
+ std::cout << "Output: " << argv[3] << std::endl;
+ std::cout << std::endl;
+
+ RLCSA rlcsa(argv[1]);
+ if(!rlcsa.supportsDisplay())
+ {
+ std::cerr << "Error: Display is not supported!" << std::endl;
+ return 2;
+ }
+ rlcsa.printInfo();
+ rlcsa.reportSize(true);
+ if(sequence >= rlcsa.getNumberOfSequences())
+ {
+ std::cerr << "Error: Invalid sequence number!" << std::endl;
+ return 3;
+ }
+
+ std::ofstream output(argv[3], std::ios_base::binary);
+ if(!output)
+ {
+ std::cerr << "Error: Cannot open output file!" << std::endl;
+ return 4;
+ }
+
+ double start = readTimer();
+ uchar* buffer = rlcsa.display(sequence);
+ usint bytes = length(rlcsa.getSequenceRange(sequence));
+ output.write((char*)buffer, bytes);
+ delete[] buffer;
+ output.close();
+
+ double time = readTimer() - start;
+ double megabytes = bytes / (double)MEGABYTE;
+ std::cout << megabytes << " megabytes in " << time << " seconds (" << (megabytes / time) << " MB/s)" << std::endl;
+ std::cout << std::endl;
+
+ return 0;
+}
--- /dev/null
+#include <algorithm>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+
+using namespace CSA;
+
+
+int main(int argc, char** argv)
+{
+ std::cout << "RLCSA LCP test" << std::endl;
+ if(argc < 5)
+ {
+ std::cout << "Usage: lcp_test basename queries runs modes [random_seed]" << std::endl;
+ std::cout << std::endl;
+ std::cout << "Supported modes:" << std::endl;
+ std::cout << "d -- Direct LCP" << std::endl;
+ std::cout << "p -- PLCP" << std::endl;
+ std::cout << "s -- Sampled LCP" << std::endl;
+ std::cout << "h -- Hybrid: PLCP and Sampled LCP" << std::endl;
+ std::cout << "l -- Locate" << std::endl;
+ std::cout << "v -- Verify results" << std::endl;
+ std::cout << std::endl;
+ return 1;
+ }
+
+ std::string base_name = argv[1];
+ std::cout << "Base name: " << base_name << std::endl;
+ usint queries = atoi(argv[2]);
+ std::cout << "Number of queries: " << queries << std::endl;
+ usint runs = std::max(1, atoi(argv[3]));
+ std::cout << "Number of test runs: " << runs << std::endl;
+
+ std::string mode_string(argv[4]);
+ usint modes = 0;
+ bool mode_direct = (mode_string.find('d') != std::string::npos);
+ bool mode_plcp = (mode_string.find('p') != std::string::npos);
+ bool mode_sampled = (mode_string.find('s') != std::string::npos);
+ bool mode_hybrid = (mode_string.find('h') != std::string::npos);
+ bool mode_locate = (mode_string.find('l') != std::string::npos);
+ bool mode_verify = (mode_string.find('v') != std::string::npos);
+ std::cout << "Modes: ";
+ if(mode_direct) { std::cout << "direct "; modes++; }
+ if(mode_plcp) { std::cout << "plcp "; modes++; }
+ if(mode_sampled) { std::cout << "sampled "; modes++; }
+ if(mode_hybrid) { std::cout << "hybrid "; modes++; }
+ if(mode_locate) { std::cout << "locate "; }
+ if(mode_verify) { std::cout << "verify"; }
+ std::cout << std::endl;
+
+ usint seed = 0xDEADBEEF;
+ if(argc > 5)
+ {
+ seed = atoi(argv[5]);
+ }
+ std::cout << "Random seed: " << seed << std::endl;
+ std::cout << std::endl;
+ if(modes == 0 || queries == 0) { return 0; }
+
+ RLCSA rlcsa(base_name);
+ if((mode_plcp || mode_hybrid) && !rlcsa.supportsLocate())
+ {
+ std::cerr << "Error: Locate is not supported!" << std::endl;
+ return 2;
+ }
+ rlcsa.printInfo();
+ rlcsa.reportSize(true);
+
+ RLEVector* plcp = 0;
+ if(mode_plcp || mode_hybrid)
+ {
+ std::string plcp_name = base_name + PLCP_EXTENSION;
+ std::ifstream plcp_file(plcp_name.c_str(), std::ios_base::binary);
+ if(!plcp_file)
+ {
+ std::cerr << "Error: Cannot open PLCP file!" << std::endl;
+ return 3;
+ }
+ plcp = new RLEVector(plcp_file);
+ std::cout << "PLCP: " << (plcp->reportSize() / (double)MEGABYTE) << " MB" << std::endl;
+ plcp_file.close();
+ }
+
+ LCPSamples* lcp = 0;
+ if(mode_sampled)
+ {
+ std::string lcp_name = base_name + LCP_SAMPLES_EXTENSION;
+ std::ifstream lcp_file(lcp_name.c_str(), std::ios_base::binary);
+ if(!lcp_file)
+ {
+ std::cerr << "Error: Cannot open LCP sample file!" << std::endl;
+ delete plcp;
+ return 4;
+ }
+ lcp = new LCPSamples(lcp_file);
+ std::cout << "Sampled LCP: " << (lcp->reportSize() / (double)MEGABYTE) << " MB" << std::endl;
+ lcp_file.close();
+ }
+
+ LCPSamples* minimal = 0;
+ if(mode_hybrid)
+ {
+ std::string minimal_name = base_name + ".minimal";
+ std::ifstream minimal_file(minimal_name.c_str(), std::ios_base::binary);
+ if(!minimal_file)
+ {
+ std::cerr << "Error: Cannot open minimal LCP sample file!" << std::endl;
+ delete plcp; delete lcp;
+ return 5;
+ }
+ minimal = new LCPSamples(minimal_file);
+ std::cout << "Minimal samples: " << (minimal->reportSize() / (double)MEGABYTE) << " MB" << std::endl;
+ minimal_file.close();
+ }
+ std::cout << std::endl;
+
+ srand(seed);
+ usint* positions = new usint[queries];
+ usint* results1 = new usint[queries];
+ usint* results2 = new usint[queries];
+ usint* results3 = new usint[queries];
+ usint* results4 = new usint[queries];
+ usint* results5 = new usint[queries];
+ for(usint i = 0; i < queries; i++)
+ {
+ positions[i] = rand() % rlcsa.getSize();
+ }
+ double start, time;
+
+ if(mode_direct)
+ {
+ std::cout << "Direct LCP computation:" << std::endl;
+ for(usint j = 0; j < runs; j++)
+ {
+ start = readTimer();
+ for(usint i = 0; i < queries; i++)
+ {
+ results1[i] = rlcsa.lcpDirect(positions[i]);
+ }
+ time = readTimer() - start;
+ std::cout << queries << " queries in " << time << " seconds (" << (queries / time) << " / s)" << std::endl;
+ }
+ std::cout << std::endl;
+ }
+
+ if(mode_plcp)
+ {
+ std::cout << "Using PLCP:" << std::endl;
+ RLEVector::Iterator iter(*plcp);
+ for(usint j = 0; j < runs; j++)
+ {
+ start = readTimer();
+ for(usint i = 0; i < queries; i++)
+ {
+ usint pos = rlcsa.locate(positions[i]);
+ results2[i] = iter.select(pos) - 2 * pos;
+ }
+ time = readTimer() - start;
+ std::cout << queries << " queries in " << time << " seconds (" << (queries / time) << " / s)" << std::endl;
+ }
+ std::cout << std::endl;
+ }
+
+ if(mode_sampled)
+ {
+ std::cout << "Using Sampled LCP:" << std::endl;
+ for(usint j = 0; j < runs; j++)
+ {
+ start = readTimer();
+ for(usint i = 0; i < queries; i++)
+ {
+ results3[i] = rlcsa.lcp(positions[i], *lcp);
+ }
+ time = readTimer() - start;
+ std::cout << queries << " queries in " << time << " seconds (" << (queries / time) << " / s)" << std::endl;
+ }
+ std::cout << std::endl;
+ }
+
+ if(mode_hybrid)
+ {
+ std::cout << "Using hybrid approach:" << std::endl;
+ for(usint j = 0; j < runs; j++)
+ {
+ start = readTimer();
+ for(usint i = 0; i < queries; i++)
+ {
+ results4[i] = rlcsa.lcp(positions[i], *minimal, *plcp);
+ }
+ time = readTimer() - start;
+ std::cout << queries << " queries in " << time << " seconds (" << (queries / time) << " / s)" << std::endl;
+ }
+ std::cout << std::endl;
+ }
+
+ if(mode_locate)
+ {
+ std::cout << "Locate:" << std::endl;
+ for(usint j = 0; j < runs; j++)
+ {
+ start = readTimer();
+ for(usint i = 0; i < queries; i++)
+ {
+ results5[i] = rlcsa.locate(positions[i]);
+ }
+ time = readTimer() - start;
+ std::cout << queries << " queries in " << time << " seconds (" << (queries / time) << " / s)" << std::endl;
+ }
+ std::cout << std::endl;
+ }
+
+ if(mode_verify && modes > 1)
+ {
+ for(usint i = 0; i < queries; i++)
+ {
+ bool ok = true;
+ ok &= !mode_direct | !mode_plcp | (results1[i] == results2[i]);
+ ok &= !mode_direct | !mode_sampled | (results1[i] == results3[i]);
+ ok &= !mode_direct | !mode_hybrid | (results1[i] == results4[i]);
+ ok &= !mode_plcp | !mode_sampled | (results2[i] == results3[i]);
+ ok &= !mode_plcp | !mode_hybrid | (results2[i] == results4[i]);
+ ok &= !mode_sampled | !mode_hybrid | (results3[i] == results4[i]);
+ if(!ok)
+ {
+ std::cout << "Query " << i << ": LCP[" << positions[i] << "] = ";
+ if(mode_direct) { std::cout << results1[i] << " (direct) "; }
+ if(mode_plcp) { std::cout << results2[i] << " (plcp) "; }
+ if(mode_sampled) { std::cout << results3[i] << " (sampled) "; }
+ if(mode_hybrid) { std::cout << results4[i] << " (hybrid) "; }
+ std::cout << std::endl;
+ }
+ }
+ }
+
+ delete[] positions;
+ delete[] results1; delete[] results2; delete[] results3; delete[] results4; delete[] results5;
+ delete plcp; delete lcp; delete minimal;
+ return 0;
+}
--- /dev/null
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+
+#include "lcpsamples.h"
+#include "misc/utils.h"
+
+
+namespace CSA
+{
+
+
+LCPSamples::LCPSamples(std::ifstream& sample_file)
+{
+ this->indexes = new DeltaVector(sample_file);
+ this->values = new Array(sample_file);
+
+ this->size = indexes->getSize();
+ this->items = indexes->getNumberOfItems();
+}
+
+LCPSamples::LCPSamples(pair_type* input, usint _size, usint _items, bool report, bool free_input) :
+ size(_size),
+ items(_items)
+{
+ DeltaEncoder index_encoder(LCPSamples::INDEX_BLOCK_SIZE);
+ ArrayEncoder value_encoder(LCPSamples::VALUE_BLOCK_SIZE);
+
+ usint max_sample = 0;
+ for(usint i = 0; i < this->items; i++)
+ {
+ index_encoder.setBit(input[i].first);
+ value_encoder.addItem(input[i].second + 1);
+ max_sample = std::max(max_sample, input[i].second);
+ }
+
+ this->indexes = new DeltaVector(index_encoder, this->size);
+ this->values = new Array(value_encoder);
+
+ if(free_input) { delete[] input; }
+
+ if(report)
+ {
+ usint max_bits = length(max_sample);
+ double total_size = max_bits * (double)(this->items) / (CHAR_BIT * (double)MEGABYTE);
+ std::cout << "Maximum sample value: " << max_sample << " (" << max_bits << " bits)" << std::endl;
+ std::cout << "Raw samples: " << total_size << " MB" << std::endl;
+ std::cout << "Encoded samples: " << (this->values->reportSize() / (double)MEGABYTE) << " MB" << std::endl;
+ std::cout << std::endl;
+ }
+}
+
+LCPSamples::LCPSamples(DeltaEncoder& index_encoder, ArrayEncoder& value_encoder, usint _size) :
+ size(_size)
+{
+ this->indexes = new DeltaVector(index_encoder, this->size);
+ this->values = new Array(value_encoder);
+ this->items = this->indexes->getNumberOfItems();
+}
+
+LCPSamples::~LCPSamples()
+{
+ delete this->indexes;
+ delete this->values;
+}
+
+void
+LCPSamples::writeTo(std::ofstream& sample_file) const
+{
+ this->indexes->writeTo(sample_file);
+ this->values->writeTo(sample_file);
+}
+
+//--------------------------------------------------------------------------
+
+usint
+LCPSamples::reportSize() const
+{
+ usint bytes = sizeof(*this);
+ bytes += this->indexes->reportSize();
+ bytes += this->values->reportSize();
+ return bytes;
+}
+
+//--------------------------------------------------------------------------
+
+
+} // namespace CSA
--- /dev/null
+#ifndef LCPSAMPLES_H
+#define LCPSAMPLES_H
+
+#include <fstream>
+
+
+#include "bits/deltavector.h"
+#include "bits/rlevector.h"
+#include "bits/nibblevector.h"
+
+#include "bits/array.h"
+#include "misc/utils.h"
+
+
+namespace CSA
+{
+
+
+class LCPSamples
+{
+ public:
+ const static usint INDEX_BLOCK_SIZE = 16;
+ const static usint VALUE_BLOCK_SIZE = 32;
+
+ LCPSamples(std::ifstream& sample_file);
+
+ // Input pairs are of form (i, LCP[i]).
+ LCPSamples(pair_type* input, usint _size, usint _items, bool report = false, bool free_input = false);
+
+ LCPSamples(DeltaEncoder& index_encoder, ArrayEncoder& value_encoder, usint _size);
+
+ ~LCPSamples();
+
+ void writeTo(std::ofstream& sample_file) const;
+
+ // Returns the value of the ith sample in suffix array order.
+ // If the index is invalid, returns this->size.
+ inline usint getSample(usint i) const
+ {
+ Array::Iterator iter(*(this->values));
+ usint tmp = iter.getItem(i);
+ if(tmp == 0) { return this->size; }
+ return tmp - 1;
+ }
+
+ // Returns (ind, sample number) where ind >= index or (size, ???).
+ inline pair_type getFirstSampleAfter(usint index) const
+ {
+ DeltaVector::Iterator iter(*(this->indexes));
+ return iter.valueAfter(index);
+ }
+
+ // Behavior is undefined if there is no sample at index.
+ inline usint getSampleAtPosition(usint index) const
+ {
+ DeltaVector::Iterator index_iter(*(this->indexes));
+ Array::Iterator value_iter(*(this->values));
+ return value_iter.getItem(index_iter.rank(index) - 1) - 1;
+ }
+
+ inline usint getNumberOfSamples() const { return this->items; }
+
+ usint reportSize() const;
+
+ private:
+ usint size, items;
+
+ DeltaVector* indexes;
+ Array* values;
+
+ // These are not allowed.
+ LCPSamples();
+ LCPSamples(const LCPSamples&);
+ LCPSamples& operator = (const LCPSamples&);
+};
+
+
+} // namespace CSA
+
+
+#endif // LCPSAMPLES_H
--- /dev/null
+#include <algorithm>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+
+using namespace CSA;
+
+
+int main(int argc, char** argv)
+{
+ std::cout << "RLCSA locate test" << std::endl;
+ if(argc < 3)
+ {
+ std::cout << "Usage: locate_test basename [begin end] output [direct]" << std::endl;
+ return 1;
+ }
+
+ usint begin = 0, end = 0;
+ int output_arg = 2;
+ std::cout << "Base name: " << argv[1] << std::endl;
+ if(argc >= 5)
+ {
+ output_arg = 4;
+ begin = atoi(argv[2]), end = atoi(argv[3]);
+ std::cout << "Begin: " << begin << std::endl;
+ std::cout << "End: " << end << std::endl;
+ if(begin > end)
+ {
+ std::cerr << "Error: Empty range!" << std::endl;
+ return 2;
+ }
+ }
+ if(argc < output_arg + 2)
+ {
+ std::cout << "Using run-based optimizations." << std::endl;
+ }
+ else
+ {
+ std::cout << "Using direct locate." << std::endl;
+ }
+ std::cout << std::endl;
+
+ RLCSA rlcsa(argv[1]);
+ if(!rlcsa.supportsLocate())
+ {
+ std::cerr << "Error: Locate is not supported!" << std::endl;
+ return 3;
+ }
+
+ if(argc >= 5)
+ {
+ if(end >= rlcsa.getSize())
+ {
+ std::cerr << "Error: Invalid range!" << std::endl;
+ return 4;
+ }
+ }
+ else
+ {
+ begin = 0; end = rlcsa.getSize() - 1;
+ }
+
+ std::ofstream output(argv[output_arg], std::ios_base::binary);
+ if(!output)
+ {
+ std::cerr << "Error: Cannot open output file!" << std::endl;
+ return 5;
+ }
+
+ usint* buffer = new usint[MILLION];
+ std::clock_t start = std::clock();
+ if(argc < output_arg + 2) // Use run-based optimizations for locate.
+ {
+ for(usint curr = begin; curr <= end; curr += MILLION)
+ {
+ pair_type range(curr, std::min(end, curr + MILLION - 1));
+ rlcsa.locate(range, buffer);
+ for(usint i = 0; i < range.second + 1 - range.first; i++)
+ {
+ output.write((char*)&(buffer[i]), sizeof(usint));
+ }
+ }
+ }
+ else // Use direct locate.
+ {
+ for(usint curr = begin; curr <= end; curr += MILLION)
+ {
+ pair_type range(curr, std::min(end, curr + MILLION - 1));
+ for(usint i = 0; i < range.second + 1 - range.first; i++)
+ {
+ buffer[i] = rlcsa.locate(curr + i);
+ }
+ for(usint i = 0; i < range.second + 1 - range.first; i++)
+ {
+ output.write((char*)&(buffer[i]), sizeof(usint));
+ }
+ }
+ }
+ std::clock_t stop = std::clock();
+ delete[] buffer;
+
+ double size = (end + 1 - begin);
+ double time = (stop - start) / (double)CLOCKS_PER_SEC;
+ std::cout << size << " locates in " << time << " seconds (" << (size / time) << " locates/s)" << std::endl;
+ output.close();
+
+ return 0;
+}
#include <climits>
+// #define MULTITHREAD_SUPPORT // Try to parallelize things using OpenMP.
+// #define MASSIVE_DATA_RLCSA // usint and sint become 64-bit in a 64-bit environment.
+
+
namespace CSA
{
typedef signed int sint;
#endif
+#ifndef uint
+typedef unsigned int uint;
+#endif
#ifndef uchar
-typedef unsigned char uchar;
+typedef unsigned char uchar;
#endif
+
typedef std::pair<usint, usint> pair_type;
return data.second + 1 - data.first;
}
+inline usint nextMultipleOf(usint multiplier, usint value)
+{
+ return multiplier * ((value / multiplier) + 1);
+}
+
const usint CHARS = ((usint)1 << CHAR_BIT);
const usint MEGABYTE = 1048576;
#include "utils.h"
+#ifdef MULTITHREAD_SUPPORT
+#include <omp.h>
+#else
+#include <cstdlib>
+#endif
+
namespace CSA
{
+//--------------------------------------------------------------------------
+
+Triple::Triple() :
+ first(0), second(0), third(0)
+{
+}
+
+Triple::Triple(usint a, usint b, usint c) :
+ first(a), second(b), third(c)
+{
+}
+
+//--------------------------------------------------------------------------
std::streamoff
fileSize(std::ifstream& file)
}
void
-readRows(std::ifstream& file, std::list<std::string>& rows, bool skipEmptyRows)
+readRows(std::ifstream& file, std::vector<std::string>& rows, bool skipEmptyRows)
{
while(file)
{
}
}
+double
+readTimer()
+{
+ #ifdef MULTITHREAD_SUPPORT
+ return omp_get_wtime();
+ #else
+ return clock() / (double)CLOCKS_PER_SEC;
+ #endif
+}
+
} // namespace CSA
#define UTILS_H
#include <fstream>
-#include <list>
+#include <vector>
#include "definitions.h"
{
+
+struct Triple
+{
+ usint first;
+ usint second;
+ usint third;
+
+ Triple();
+ Triple(usint a, usint b, usint c);
+};
+
+
std::streamoff fileSize(std::ifstream& file);
std::streamoff fileSize(std::ofstream& file);
std::ostream& operator<<(std::ostream& stream, pair_type data);
-void readRows(std::ifstream& file, std::list<std::string>& rows, bool skipEmptyRows);
+void readRows(std::ifstream& file, std::vector<std::string>& rows, bool skipEmptyRows);
+
+double readTimer();
} // namespace CSA
--- /dev/null
+#include <algorithm>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+#ifdef MULTITHREAD_SUPPORT
+#include <omp.h>
+#endif
+
+#include "rlcsa_builder.h"
+#include "misc/utils.h"
+
+
+using namespace CSA;
+
+
+double indexParts(std::vector<std::string>& filename, usint threads, Parameters& parameters);
+
+const int MAX_THREADS = 64;
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "Parallel RLCSA builder" << std::endl;
+ if(argc < 3)
+ {
+ std::cout << "Usage: parallel_build [-n] listname output [threads]" << std::endl;
+ std::cout << " -n do not merge the indexes" << std::endl;
+ return 1;
+ }
+
+ int list_parameter = 1, output_parameter = 2, threads_parameter = 3;
+ bool do_merge = true;
+ if(std::string("-n").compare(argv[1]) == 0)
+ {
+ list_parameter++; output_parameter++; threads_parameter++;
+ do_merge = false;
+ std::cout << "Option '-n' specified. Partial indexes will not be merged." << std::endl;
+ }
+
+ std::ifstream filelist(argv[list_parameter], std::ios_base::binary);
+ if(!filelist)
+ {
+ std::cerr << "Error opening file list!" << std::endl;
+ return 2;
+ }
+ std::vector<std::string> files;
+ readRows(filelist, files, true);
+ filelist.close();
+ std::cout << "Input files: " << files.size() << std::endl;
+
+ std::string base_name = argv[output_parameter];
+ std::cout << "Output: " << base_name << std::endl;
+
+ usint threads = 1;
+ if(argc > threads_parameter)
+ {
+ threads = std::min(MAX_THREADS, std::max(atoi(argv[threads_parameter]), 1));
+ }
+ std::cout << "Threads: " << threads << std::endl;
+ std::cout << std::endl;
+
+ std::string parameters_name = base_name + PARAMETERS_EXTENSION;
+ Parameters parameters;
+ parameters.set(RLCSA_BLOCK_SIZE);
+ parameters.set(SAMPLE_RATE);
+ parameters.set(SUPPORT_LOCATE);
+ parameters.set(SUPPORT_DISPLAY);
+ parameters.read(parameters_name);
+ parameters.print();
+
+ double start = readTimer();
+ double megabytes = indexParts(files, threads, parameters);
+
+ RLCSABuilder builder(parameters.get(RLCSA_BLOCK_SIZE), parameters.get(SAMPLE_RATE), 0, threads);
+ if(do_merge)
+ {
+ std::cout << "Phase 2: Merging the indexes" << std::endl;
+ for(std::vector<std::string>::iterator iter = files.begin(); iter != files.end(); iter++)
+ {
+ std::cout << "Increment: " << *iter << std::endl;
+ builder.insertFromFile(*iter);
+ }
+ std::cout << std::endl;
+
+ RLCSA* index = builder.getRLCSA();
+ if(index != 0 && index->isOk())
+ {
+ index->printInfo();
+ index->reportSize(true);
+ index->writeTo(base_name);
+ parameters.write(parameters_name);
+ }
+ delete index;
+ }
+
+ double stop = readTimer();
+ std::cout << megabytes << " megabytes indexed in " << (stop - start) << " seconds (" << (megabytes / (stop - start)) << " MB/s)." << std::endl;
+ if(do_merge)
+ {
+ std::cout << "Search time: " << builder.getSearchTime() << " seconds" << std::endl;
+ std::cout << "Sort time: " << builder.getSortTime() << " seconds" << std::endl;
+ std::cout << "Merge time: " << builder.getMergeTime() << " seconds" << std::endl;
+ }
+ std::cout << std::endl;
+
+ return 0;
+}
+
+
+double
+indexParts(std::vector<std::string>& filenames, usint threads, Parameters& parameters)
+{
+ double start = readTimer();
+ std::cout << "Phase 1: Building indexes for input files" << std::endl;
+ usint block_size = parameters.get(RLCSA_BLOCK_SIZE);
+ usint sample_rate = parameters.get(SAMPLE_RATE);
+ usint total_size = 0;
+
+ std::ifstream* input_file;
+ usint size;
+ std::string parameters_name;
+ RLCSA* index;
+ uchar* data;
+ sint i;
+
+ #ifdef MULTITHREAD_SUPPORT
+ omp_set_num_threads(threads);
+ #pragma omp parallel private(input_file, size, parameters_name, index, data)
+ {
+ #pragma omp for schedule(dynamic, 1)
+ #endif
+ for(i = 0; i < (sint)(filenames.size()); i++)
+ {
+ #ifdef MULTITHREAD_SUPPORT
+ #pragma omp critical
+ {
+ #endif
+ size = 0; data = 0;
+ std::cout << "Input: " << filenames[i] << std::endl;
+ input_file = new std::ifstream(filenames[i].c_str(), std::ios_base::binary);
+ if(input_file == 0)
+ {
+ std::cerr << "Error opening input file " << filenames[i] << "!" << std::endl;
+ }
+ else
+ {
+ size = fileSize(*input_file);
+ data = new uchar[size];
+ input_file->read((char*)data, size);
+ delete input_file;
+ }
+ #ifdef MULTITHREAD_SUPPORT
+ }
+ #endif
+
+ if(size > 0)
+ {
+ index = new RLCSA(data, size, block_size, sample_rate, true, true);
+ if(index != 0 && index->isOk()) { index->writeTo(filenames[i]); }
+ delete index;
+
+ #ifdef MULTITHREAD_SUPPORT
+ #pragma omp critical
+ {
+ #endif
+ total_size += size;
+ parameters_name = filenames[i] + PARAMETERS_EXTENSION;
+ parameters.write(parameters_name);
+ #ifdef MULTITHREAD_SUPPORT
+ }
+ #endif
+ }
+ else
+ {
+ #ifdef MULTITHREAD_SUPPORT
+ #pragma omp critical
+ {
+ #endif
+ std::cerr << "Warning: Empty input file " << filenames[i] << "!" << std::endl;
+ #ifdef MULTITHREAD_SUPPORT
+ }
+ #endif
+ }
+ }
+ #ifdef MULTITHREAD_SUPPORT
+ }
+ #endif
+
+ double stop = readTimer();
+ double megabytes = total_size / (double)MEGABYTE;
+ std::cout << "Indexed " << megabytes << " megabytes in " << (stop - start) << " seconds (" << (megabytes / (stop - start)) << " MB/s)." << std::endl;
+ std::cout << std::endl;
+
+ return megabytes;
+}
of this software.*/
/*
- Replaced int -> sint for massive data support and moved into namespace CSA.
+ Moved everything into a struct to support parallel construction.
+ Created another version using sint instead of int to index more than 2 GB at a time.
+ Put everything into a namespace.
- - Jouni Siren 2009-03-16
+ - Jouni Siren 2009-03-20
*/
#include <climits>
{
-static sint *I, /* group array, ultimately suffix array.*/
+#define KEY(p) (V[*(p)+(h)])
+#define SWAP(p, q) (tmp=*(p), *(p)=*(q), *(q)=tmp)
+#define MED3(a, b, c) (KEY(a)<KEY(b) ? (KEY(b)<KEY(c) ? b : KEY(a)<KEY(c) ? c : (a)) : (KEY(b)>KEY(c) ? b : KEY(a)>KEY(c) ? c : (a)))
+
+
+//--------------------------------------------------------------------------
+
+template <class T>
+struct QSufSort
+{
+
+ T *I, /* group array, ultimately suffix array.*/
*V, /* inverse array, ultimately inverse of I.*/
r, /* number of symbols aggregated by transform.*/
h; /* length of already-sorted prefixes.*/
-#define KEY(p) (V[*(p)+(h)])
-#define SWAP(p, q) (tmp=*(p), *(p)=*(q), *(q)=tmp)
-#define MED3(a, b, c) (KEY(a)<KEY(b) ? (KEY(b)<KEY(c) ? b : KEY(a)<KEY(c) ? c : (a)) : (KEY(b)>KEY(c) ? b : KEY(a)>KEY(c) ? c : (a)))
/* Subroutine for select_sort_split and sort_split. Sets group numbers for a
group whose lowest position in I is pl and highest position is pm.*/
-static void update_group(sint *pl, sint *pm)
+void update_group(T *pl, T *pm)
{
- sint g;
+ T g;
g=pm-I; /* group number.*/
V[*pl]=g; /* update group number of first position.*/
/* Quadratic sorting method to use for small subarrays. To be able to update
group numbers consistently, a variant of selection sorting is used.*/
-static void select_sort_split(sint *p, sint n) {
- sint *pa, *pb, *pi, *pn;
- sint f, v, tmp;
+void select_sort_split(T *p, T n) {
+ T *pa, *pb, *pi, *pn;
+ T f, v, tmp;
pa=p; /* pa is start of group being picked out.*/
pn=p+n-1; /* pn is last position of subarray.*/
/* Subroutine for sort_split, algorithm by Bentley & McIlroy.*/
-static sint choose_pivot(sint *p, sint n) {
- sint *pl, *pm, *pn;
- sint s;
+T choose_pivot(T *p, T n) {
+ T *pl, *pm, *pn;
+ T s;
pm=p+(n>>1); /* small arrays, middle element.*/
if (n>7) {
Software -- Practice and Experience 23(11), 1249-1265 (November 1993). This
function is based on Program 7.*/
-static void sort_split(sint *p, sint n)
+void sort_split(T *p, T n)
{
- sint *pa, *pb, *pc, *pd, *pl, *pm, *pn;
- sint f, v, s, t, tmp;
+ T *pa, *pb, *pc, *pd, *pl, *pm, *pn;
+ T f, v, s, t, tmp;
if (n<7) { /* multi-selection sort smallest arrays.*/
select_sort_split(p, n);
Output: x is V and p is I after the initial sorting stage of the refined
suffix sorting algorithm.*/
-
-static void bucketsort(sint *x, sint *p, sint n, sint k)
+
+void bucketsort(T *x, T *p, T n, T k)
{
- sint *pi, i, c, d, g;
+ T *pi, i, c, d, g;
for (pi=p; pi<p+k; ++pi)
*pi=-1; /* mark linked lists empty.*/
new alphabet. If j<=n+1, the alphabet is compacted. The global variable r is
set to the number of old symbols grouped into one. Only x[n] is 0.*/
-static sint transform(sint *x, sint *p, sint n, sint k, sint l, sint q)
+T transform(T *x, T *p, T n, T k, T l, T q)
{
- sint b, c, d, e, i, j, m, s;
- sint *pi, *pj;
+ T b, c, d, e, i, j, m, s;
+ T *pi, *pj;
for (s=0, i=k-l; i; i>>=1)
++s; /* s is number of bits in old symbol.*/
contents of x[n] is disregarded, the n-th symbol being regarded as
end-of-string smaller than all other symbols.*/
-void suffixsort(sint *x, sint *p, sint n, sint k, sint l)
+void suffixsort(T *x, T *p, T n, T k, T l)
{
- sint *pi, *pk;
- sint i, j, s, sl;
+ T *pi, *pk;
+ T i, j, s, sl;
V=x; /* set global values.*/
I=p;
I[V[i]]=i;
}
+}; // QSufSort
+
+//--------------------------------------------------------------------------
+
+
+void suffixsort(int *x, int *p, int n, int k, int l)
+{
+ QSufSort<int> sorter;
+ sorter.suffixsort(x, p, n, k, l);
+}
+
+void massive_suffixsort(sint *x, sint *p, sint n, sint k, sint l)
+{
+ QSufSort<sint> sorter;
+ sorter.suffixsort(x, p, n, k, l);
+}
+
} // namespace CSA
{
-void suffixsort(sint *, sint *, sint, sint, sint);
+void suffixsort(int *, int *, int, int, int);
+
+void massive_suffixsort(sint *, sint *, sint, sint, sint);
} // namespace CSA
--- /dev/null
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+
+
+using namespace CSA;
+
+
+/*
+ This program writes run-length encoded PLCP of the collection into a file.
+*/
+
+
+usint
+countRuns(usint prev, uchar* buffer, usint length)
+{
+ usint runs = 0;
+
+ for(usint i = 0; i < length; i++)
+ {
+ if(buffer[i] != prev)
+ {
+ prev = buffer[i];
+ if(buffer[i] != 0) { runs++; }
+ }
+ }
+
+ return runs;
+}
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "RLCSA to BWT converter" << std::endl;
+ if(argc < 2)
+ {
+ std::cout << "Usage: read_bwt base_name [buffer_size]" << std::endl;
+ return 1;
+ }
+
+ std::string base_name = argv[1];
+ std::string bwt_name = base_name + ".bwt";
+ std::cout << "BWT: " << bwt_name << std::endl;
+ std::ofstream bwt_file(bwt_name.c_str(), std::ios_base::binary);
+ if(!bwt_file)
+ {
+ std::cerr << "Error creating BWT file!" << std::endl;
+ return 2;
+ }
+ std::cout << std::endl;
+
+ RLCSA rlcsa(base_name);
+ clock_t start = clock();
+ usint buffer_size = 0;
+ if(argc > 2) { buffer_size = atoi(argv[2]); }
+ usint n = rlcsa.getSize() + rlcsa.getNumberOfSequences();
+
+ usint runs = 0, prev = CHARS;
+ if(buffer_size > 0)
+ {
+ for(usint i = 0; i < n; i += buffer_size)
+ {
+ pair_type range(i, std::min(i + buffer_size - 1, n - 1));
+ uchar* bwt = rlcsa.readBWT(range);
+ if(bwt != 0)
+ {
+ runs += countRuns(prev, bwt, length(range));
+ prev = bwt[length(range) - 1];
+ bwt_file.write((char*)bwt, length(range));
+ delete[] bwt;
+ }
+ }
+ }
+ else
+ {
+ uchar* bwt = rlcsa.readBWT();
+ if(bwt != 0)
+ {
+ runs = countRuns(prev, bwt, n);
+ bwt_file.write((char*)bwt, n);
+ delete[] bwt;
+ }
+ }
+
+ clock_t stop = clock();
+ double time = ((stop - start) / (double)CLOCKS_PER_SEC);
+ double megabytes = n / (double)MEGABYTE;
+ std::cout << megabytes << " megabytes in " << time << " seconds (" << (megabytes / time) << " MB/s)" << std::endl;
+ std::cout << std::endl;
+
+ // Testing direct reporting of the number of runs.
+ // This is as expensive as reading the BWT.
+ std::cout << "Number of runs: " << runs << std::endl;
+ runs = rlcsa.countRuns();
+ std::cout << "Number of runs (direct count): " << runs << std::endl;
+ std::cout << std::endl;
+
+ bwt_file.close();
+ return 0;
+}
#include <algorithm>
#include <cstdlib>
-#include <ctime>
+#include <cstring>
#include <fstream>
#include <iostream>
+#include <list>
#include "rlcsa.h"
#include "misc/utils.h"
#include "qsufsort/qsufsort.h"
+#include "bits/vectors.h"
+
+
+#ifdef MULTITHREAD_SUPPORT
+#include <omp.h>
+#endif
+
namespace CSA
array_file.read((char*)distribution, CHARS * sizeof(usint));
this->buildCharIndexes(distribution);
- Parameters parameters;
- parameters.read(base_name + PARAMETERS_EXTENSION);
+// Parameters parameters;
+// parameters.read(base_name + PARAMETERS_EXTENSION);
for(usint c = 0; c < CHARS; c++)
{
- if(!isEmpty(this->index_ranges[c])) { this->array[c] = new RLEVector(array_file); }
+ if(!isEmpty(this->index_ranges[c])) { this->array[c] = new PsiVector(array_file); }
}
this->end_points = new DeltaVector(array_file);
array_file.read((char*)&(this->sample_rate), sizeof(this->sample_rate));
array_file.close();
- this->support_locate = parameters.get(SUPPORT_LOCATE);
- this->support_display = parameters.get(SUPPORT_DISPLAY);
+ this->support_locate = true;// parameters.get(SUPPORT_LOCATE);
+ this->support_display = true; //parameters.get(SUPPORT_DISPLAY);
if(this->support_locate || this->support_display)
{
sa_sample_file.close();
}
- if(print) { parameters.print(); }
+// if(print) { parameters.print(); }
this->ok = true;
}
// Build suffix array.
- sint* inverse = new sint[bytes + 1];
+ int* inverse = new int[bytes + 1];
if(multiple_sequences)
{
- sint zeros = 0;
+ int zeros = 0;
for(usint i = 0; i < bytes; i++)
{
if(data[i] == 0)
}
else
{
- inverse[i] = (sint)data[i] + this->number_of_sequences;
+ inverse[i] = (int)(data[i] + this->number_of_sequences);
}
}
}
else
{
- for(usint i = 0; i < bytes; i++) { inverse[i] = (sint)data[i]; }
+ for(usint i = 0; i < bytes; i++) { inverse[i] = (int)data[i]; }
}
if(delete_data) { delete[] data; }
- sint* sa = new sint[bytes + 1];
+ int* sa = new int[bytes + 1];
suffixsort(inverse, sa, bytes, high + 1, low);
// Sample SA.
- usint incr = (multiple_sequences ? this->number_of_sequences + 1 : 1); // No e_of_s markers in SA.
if(this->support_locate || this->support_display)
{
if(multiple_sequences)
{
- std::cout << "We shouldn't be here!" << std::endl;
- // Real SA starts at sa + incr.
- // for each sequence
- // find starting position
- // sample relative to starting position
- // sort samples to SA order
+ this->sa_samples = new SASamples((uint*)inverse, this->end_points, this->data_size, this->sample_rate);
}
else
{
- this->sa_samples = new SASamples((usint*)(sa + incr), this->data_size, this->sample_rate);
+ this->sa_samples = new SASamples((uint*)(sa + 1), this->data_size, this->sample_rate);
}
}
// Build RLCSA.
+ usint incr = (multiple_sequences ? this->number_of_sequences + 1 : 1); // No e_of_s markers in SA.
usint decr = (multiple_sequences ? 1 : 0); // Ignore the last e_of_s marker if multiple sequences.
for(usint c = 0; c < CHARS; c++)
{
if(distribution[c] == 0) { this->array[c] = 0; continue; }
- usint* curr = (usint*)(sa + index_ranges[c].first + incr);
- usint* limit = (usint*)(sa + index_ranges[c].second + incr + 1);
- RLEEncoder encoder(block_size);
+ uint* curr = (uint*)(sa + index_ranges[c].first + incr);
+ uint* limit = (uint*)(sa + index_ranges[c].second + incr + 1);
+ PsiVector::Encoder encoder(block_size);
pair_type run(*curr, 1);
curr++;
}
encoder.setRun(run.first - decr, run.second);
- this->array[c] = new RLEVector(encoder, this->data_size + incr - decr);
+ this->array[c] = new PsiVector(encoder, this->data_size + incr - decr);
}
delete[] sa;
this->ok = true;
}
-RLCSA::RLCSA(RLCSA& index, RLCSA& increment, usint* positions, usint block_size) :
+RLCSA::RLCSA(RLCSA& index, RLCSA& increment, usint* positions, usint block_size, usint threads) :
ok(false),
sa_samples(0),
end_points(0)
if(index.sample_rate != increment.sample_rate)
{
std::cerr << "RLCSA: Cannot combine indexes with different sample rates!" << std::endl;
- std::cout << "Index: " << index.sample_rate << ", increment: " << increment.sample_rate << std::endl;
return;
}
}
this->buildCharIndexes(distribution);
this->sample_rate = index.sample_rate;
-
-
- // Combine end points of sequences.
this->number_of_sequences = index.number_of_sequences + increment.number_of_sequences;
- DeltaEncoder* endings = new DeltaEncoder(RLCSA::ENDPOINT_BLOCK_SIZE);
-
- endings->setBit(index.end_points->select(0));
- for(usint i = 1; i < index.number_of_sequences; i++)
- {
- endings->setBit(index.end_points->selectNext());
- }
- usint sum = index.end_points->getSize();
- delete index.end_points; index.end_points = 0;
-
- endings->setBit(sum + increment.end_points->select(0));
- for(usint i = 1; i < increment.number_of_sequences; i++)
- {
- endings->setBit(sum + increment.end_points->selectNext());
- }
- sum += increment.end_points->getSize();
- delete increment.end_points; increment.end_points = 0;
-
- this->end_points = new DeltaVector(*endings, sum);
- delete endings;
-
- // Combine Psi.
+ // Merge end points, SA samples, and Psi.
usint psi_size = this->data_size + this->number_of_sequences;
- for(usint c = 0; c < CHARS; c++)
- {
- if(distribution[c] == 0) { this->array[c] = 0; continue; }
- this->array[c] = mergeVectors(index.array[c], increment.array[c], positions, increment.data_size + increment.number_of_sequences, psi_size, block_size);
- index.array[c] = 0;
- increment.array[c] = 0;
-
- if(this->array[c] == 0)
- {
- std::cerr << "RLCSA: Merge failed for vectors " << c << "!" << std::endl;
- return;
- }
- }
-
+ bool should_be_ok = true;
- // Combine suffix array samples.
- if(this->support_locate || this->support_display)
+ #ifdef MULTITHREAD_SUPPORT
+ omp_set_num_threads(threads);
+ #pragma omp parallel
{
- positions += increment.number_of_sequences;
- for(usint i = 0; i < increment.data_size; i++)
+ #pragma omp for schedule(dynamic, 1)
+ #endif
+ for(int c = -2; c < (int)CHARS; c++)
{
- positions[i] -= this->number_of_sequences;
+ if(c == -2) { this->mergeEndPoints(index, increment); }
+ else if(c == -1) { this->mergeSamples(index, increment, positions); }
+ else if(distribution[c] != 0)
+ {
+ this->array[c] = mergeVectors<PsiVector, PsiVector::Encoder, PsiVector::Iterator>(index.array[c], increment.array[c], positions, increment.data_size + increment.number_of_sequences, psi_size, block_size);
+ index.array[c] = 0;
+ increment.array[c] = 0;
+
+ if(this->array[c] == 0)
+ {
+ std::cerr << "RLCSA: Merge failed for vectors " << c << "!" << std::endl;
+ should_be_ok = false;
+ }
+ }
}
- this->sa_samples = new SASamples(*(index.sa_samples), *(increment.sa_samples), positions, increment.data_size);
+ #ifdef MULTITHREAD_SUPPORT
}
+ #endif
- this->ok = true;
+ this->ok = should_be_ok;
}
RLCSA::~RLCSA()
delete this->end_points;
}
+//--------------------------------------------------------------------------
+
void
-RLCSA::writeTo(const std::string& base_name)
+RLCSA::writeTo(const std::string& base_name) const
{
std::string array_name = base_name + ARRAY_EXTENSION;
std::ofstream array_file(array_name.c_str(), std::ios_base::binary);
}
bool
-RLCSA::isOk()
+RLCSA::isOk() const
{
return this->ok;
}
//--------------------------------------------------------------------------
pair_type
-RLCSA::count(const std::string& pattern)
+RLCSA::count(const std::string& pattern) const
{
if(pattern.length() == 0) { return pair_type(0, this->data_size - 1); }
- pair_type index_range = this->index_ranges[(usint)*(pattern.rbegin())];
- if(isEmpty(index_range)) { return index_range; }
+ std::string::const_reverse_iterator iter = pattern.rbegin();
+
+ pair_type index_range = this->index_ranges[(uchar)*iter];
+ if(isEmpty(index_range)) { return EMPTY_PAIR; }
index_range.first += this->number_of_sequences;
index_range.second += this->number_of_sequences;
- for(std::string::const_reverse_iterator iter = ++pattern.rbegin(); iter != pattern.rend(); iter++)
+ for(++iter; iter != pattern.rend(); ++iter)
{
- RLEVector* vector = this->array[(usint)*iter];
- usint start = this->index_ranges[(usint)*iter].first;
-
- index_range.first = start + vector->rank(index_range.first, true) - 1 + this->number_of_sequences;
- index_range.second = start + vector->rank(index_range.second) - 1 + this->number_of_sequences;
-
- if(isEmpty(index_range)) { return index_range; }
+ index_range = this->backwardSearchStep(index_range, (uchar)*iter);
+ if(isEmpty(index_range)) { return EMPTY_PAIR; }
}
// Suffix array indexes are 0-based.
return index_range;
}
+//--------------------------------------------------------------------------
+
void
-RLCSA::reportPositions(uchar* data, usint length, usint* positions)
+RLCSA::reportPositions(uchar* data, usint length, usint* positions) const
{
if(data == 0 || length == 0 || positions == 0) { return; }
+ PsiVector::Iterator** iters = this->getIterators();
+
usint current = this->number_of_sequences - 1;
positions[length] = current; // "immediately after current"
for(sint i = (sint)(length - 1); i >= 0; i--)
{
-// positions[i] = current; // "immediately after current"
usint c = (usint)data[i];
- if(array[c])
+ if(this->array[c] != 0)
{
- current = this->index_ranges[c].first + this->array[c]->rank(current) - 1 + this->number_of_sequences;
+ current = this->LF(current, c, *(iters[c]));
}
else
{
}
positions[i] = current; // "immediately after current"
}
+
+ this->deleteIterators(iters);
}
//--------------------------------------------------------------------------
-LocateItem*
-RLCSA::locate(pair_type range)
+usint*
+RLCSA::locate(pair_type range) const
{
if(!(this->support_locate) || isEmpty(range) || range.second >= this->data_size) { return 0; }
- LocateItem* data = new LocateItem[range.second + 1 - range.first];
- if(!data) { return 0; }
+ usint* data = new usint[length(range)];
this->locateUnsafe(range, data);
return data;
}
-LocateItem*
-RLCSA::locate(pair_type range, LocateItem* data)
+usint*
+RLCSA::locate(pair_type range, usint* data) const
{
if(!(this->support_locate) || isEmpty(range) || range.second >= this->data_size || data == 0) { return 0; }
this->locateUnsafe(range, data);
}
void
-RLCSA::locateUnsafe(pair_type range, LocateItem* data)
+RLCSA::locateUnsafe(pair_type range, usint* data) const
{
- usint items = range.second + 1 - range.first;
+ usint items = length(range);
+ usint* offsets = new usint[items];
+ bool* finished = new bool[items]; // FIXME This could be more space efficient...
+
+ PsiVector::Iterator** iters = this->getIterators();
+
for(usint i = 0, j = range.first; i < items; i++, j++)
{
- data[i].value = j + this->number_of_sequences;
- data[i].offset = 0;
- data[i].found = false;
+ data[i] = j + this->number_of_sequences;
+ offsets[i] = 0;
+ finished[i] = false;
}
bool found = false;
pair_type run = EMPTY_PAIR;
for(usint i = 0; i < items; i++)
{
- if(data[i].found)
+ if(finished[i])
{
continue; // The run might continue after this.
}
{
run = pair_type(i, i);
}
- else if(data[i].value - data[run.first].value == i - run.first)
+ else if(data[i] - data[run.first] == i - run.first)
{
run.second = i;
}
else
{
- found &= this->processRun(run, data);
+ found &= this->processRun(run, data, offsets, finished, iters);
run = pair_type(i, i);
}
}
- if(!isEmpty(run)) { found &= this->processRun(run, data); }
+ if(!isEmpty(run)) { found &= this->processRun(run, data, offsets, finished, iters); }
}
+
+ this->deleteIterators(iters);
+ delete[] offsets;
+ delete[] finished;
}
bool
-RLCSA::processRun(pair_type run, LocateItem* data)
+RLCSA::processRun(pair_type run, usint* data, usint* offsets, bool* finished, PsiVector::Iterator** iters) const
{
bool found = true;
usint run_start = 0, run_left = 0;
for(usint i = run.first; i <= run.second; i++)
{
- if(data[i].found)
+ if(finished[i])
{
if(run_left > 0) { run_left--; }
continue;
}
- if(data[i].value < this->number_of_sequences) // Implicit sample here.
+ if(data[i] < this->number_of_sequences) // Implicit sample here.
{
- data[i].value = this->end_points->select(data[i].value) + 1 - data[i].offset;
- data[i].found = true;
+ DeltaVector::Iterator iter(*(this->end_points));
+ data[i] = iter.select(data[i]) + 1 - offsets[i];
+ finished[i] = true;
if(run_left > 0) { run_left--; }
continue;
}
- if(next_sample.first < data[i].value) // Need another sample.
+ if(next_sample.first < data[i]) // Need another sample.
{
- next_sample = this->sa_samples->getFirstSampleAfter(data[i].value - this->number_of_sequences);
+ next_sample = this->sa_samples->getFirstSampleAfter(data[i] - this->number_of_sequences);
next_sample.first += this->number_of_sequences;
}
- if(data[i].value < next_sample.first) // No sample found for current position.
+ if(data[i] < next_sample.first) // No sample found for current position.
{
if(run_left > 0)
{
- data[i].value = data[run_start].value + i - run_start;
+ data[i] = data[run_start] + i - run_start;
run_left--;
}
else
{
- pair_type value = this->psi(data[i].value - this->number_of_sequences, run.second - i);
- data[i].value = value.first;
+ pair_type value = this->psi(data[i] - this->number_of_sequences, run.second - i, iters);
+ data[i] = value.first;
run_left = value.second;
run_start = i;
}
- data[i].offset++;
+ offsets[i]++;
found = false;
}
else // Sampled position found.
{
- data[i].value = this->sa_samples->getSample(next_sample.second) - data[i].offset;
- data[i].found = true;
+ data[i] = this->sa_samples->getSample(next_sample.second) - offsets[i];
+ finished[i] = true;
if(run_left > 0) { run_left--; }
}
}
}
usint
-RLCSA::locate(usint index)
+RLCSA::locate(usint index) const
{
if(!(this->support_locate) || index >= this->data_size) { return this->data_size; }
{
if(index < this->number_of_sequences) // Implicit sample here
{
- return this->end_points->select(index) + 1 - offset;
+ DeltaVector::Iterator iter(*(this->end_points));
+ return iter.select(index) + 1 - offset;
}
pair_type next_sample = this->sa_samples->getFirstSampleAfter(index - this->number_of_sequences);
next_sample.first += this->number_of_sequences;
//--------------------------------------------------------------------------
uchar*
-RLCSA::display(usint sequence, pair_type range)
+RLCSA::display(usint sequence) const
+{
+ if(!(this->support_display)) { return 0; }
+
+ pair_type seq_range = this->getSequenceRange(sequence);
+ if(isEmpty(seq_range)) { return 0; }
+
+ uchar* data = new uchar[length(seq_range)+1];
+ this->displayUnsafe(seq_range, data);
+ data[length(seq_range)] = 0;
+
+ return data;
+}
+
+uchar*
+RLCSA::display(usint sequence, pair_type range) const
{
if(!(this->support_display) || isEmpty(range)) { return 0; }
- pair_type seq_range = this->getSequence(sequence);
+ pair_type seq_range = this->getSequenceRange(sequence);
if(isEmpty(seq_range)) { return 0; }
range.first += seq_range.first; range.second += seq_range.first;
if(range.second > seq_range.second) { return 0; }
- uchar* data = new uchar[range.second + 1 - range.first];
- if(!data) { return 0; }
+ uchar* data = new uchar[length(range)];
this->displayUnsafe(range, data);
return data;
}
uchar*
-RLCSA::display(usint sequence, pair_type range, uchar* data)
+RLCSA::display(usint sequence, pair_type range, uchar* data) const
{
if(!(this->support_display) || isEmpty(range) || data == 0) { return 0; }
- pair_type seq_range = this->getSequence(sequence);
+ pair_type seq_range = this->getSequenceRange(sequence);
if(isEmpty(seq_range)) { return 0; }
range.first += seq_range.first; range.second += seq_range.first;
return data;
}
+uchar*
+RLCSA::display(usint position, usint len, usint context, usint& result_length) const
+{
+ if(!(this->support_display)) { return 0; }
+
+ pair_type range = this->getSequenceRangeForPosition(position);
+ if(isEmpty(range)) { return 0; }
+
+ range.first = position - std::min(context, position - range.first);
+ range.second = std::min(range.second, position + len + context - 1);
+ result_length = length(range);
+ if(isEmpty(range)) { return 0; }
+
+ uchar* data = new uchar[length(range)];
+ this->displayUnsafe(range, data);
+
+ return data;
+}
+
+usint
+RLCSA::displayPrefix(usint sequence, usint len, uchar* data) const
+{
+ if(!(this->support_display) || len == 0 || data == 0) { return 0; }
+
+ pair_type seq_range = this->getSequenceRange(sequence);
+ if(isEmpty(seq_range)) { return 0; }
+
+ pair_type range(seq_range.first, std::min(seq_range.second, seq_range.first + len - 1));
+
+ this->displayUnsafe(range, data);
+ return length(range);
+}
+
+usint
+RLCSA::displayFromPosition(usint index, usint max_len, uchar* data) const
+{
+ if(max_len == 0 || data == 0 || index >= this->data_size) { return 0; }
+
+ for(usint i = 0; i < max_len; i++)
+ {
+ data[i] = this->getCharacter(index);
+ index = this->psiUnsafe(index, data[i]);
+ if(index < this->number_of_sequences) { return i + 1; }
+ index -= this->number_of_sequences;
+ }
+
+ return max_len;
+}
+
void
-RLCSA::displayUnsafe(pair_type range, uchar* data)
+RLCSA::displayUnsafe(pair_type range, uchar* data) const
{
usint i = range.first - range.first % this->sa_samples->getSampleRate();
usint pos = this->sa_samples->inverseSA(i);
- for(; i < range.first; i++)
+ if(length(range) >= 1024)
{
- pos = this->psi(pos) - this->number_of_sequences;
+ PsiVector::Iterator** iters = this->getIterators();
+ for(; i < range.first; i++)
+ {
+ pos = this->psi(pos, iters) - this->number_of_sequences;
+ }
+ for(; i <= range.second; i++)
+ {
+ usint c = this->getCharacter(pos);
+ data[i - range.first] = c;
+ pos = this->psiUnsafe(pos, c, *(iters[c])) - this->number_of_sequences;
+ }
+ this->deleteIterators(iters);
}
- for(; i <= range.second; i++)
+ else
{
- data[i - range.first] = this->getCharacter(pos);
- pos = this->psi(pos) - this->number_of_sequences;
+ for(; i < range.first; i++)
+ {
+ pos = this->psi(pos) - this->number_of_sequences;
+ }
+ for(; i <= range.second; i++)
+ {
+ usint c = this->getCharacter(pos);
+ data[i - range.first] = c;
+ pos = this->psiUnsafe(pos, c) - this->number_of_sequences;
+ }
}
}
//--------------------------------------------------------------------------
-void
-RLCSA::decompressInto(std::ofstream& psi_file)
+pair_type
+RLCSA::getSequenceRange(usint number) const
{
- for(usint c = 0; c < CHARS; c++)
+ if(number >= this->number_of_sequences) { return EMPTY_PAIR; }
+
+ pair_type result;
+ DeltaVector::Iterator iter(*(this->end_points));
+ if(number == 0)
{
- if(!(this->array[c])) { continue; }
- usint value = this->array[c]->select(0);
- psi_file.write((char*)&(value), sizeof(value));
- while(this->array[c]->hasNext())
- {
- value = this->array[c]->selectNext();
- psi_file.write((char*)&value, sizeof(value));
- }
+ result.first = 0;
+ result.second = iter.select(number);
+ }
+ else
+ {
+ result.first = nextMultipleOf(this->sample_rate, iter.select(number - 1));
+ result.second = iter.selectNext();
+ }
+
+ return result;
+}
+
+usint
+RLCSA::getSequenceForPosition(usint value) const
+{
+ if(value == 0) { return 0; }
+ DeltaVector::Iterator iter(*(this->end_points));
+ return iter.rank(value - 1);
+}
+
+pair_type
+RLCSA::getSequenceRangeForPosition(usint value) const
+{
+ return this->getSequenceRange(this->getSequenceForPosition(value));
+}
+
+usint*
+RLCSA::getSequenceForPosition(usint* values, usint len) const
+{
+ if(values == 0) { return 0; }
+
+ DeltaVector::Iterator iter(*(this->end_points));
+ for(usint i = 0; i < len; i++)
+ {
+ if(values[i] > 0) { values[i] = iter.rank(values[i] - 1); }
+ }
+
+ return values;
+}
+
+pair_type
+RLCSA::getRelativePosition(usint value) const
+{
+ DeltaVector::Iterator iter(*(this->end_points));
+ pair_type result(0, value);
+
+ if(value > 0) { result.first = iter.rank(value - 1); }
+ if(result.first > 0)
+ {
+ result.second -= nextMultipleOf(this->sample_rate, iter.select(result.first - 1));
}
+
+ return result;
+}
+
+//--------------------------------------------------------------------------
+
+uchar*
+RLCSA::readBWT() const
+{
+ return this->readBWT(pair_type(0, this->data_size + this->number_of_sequences - 1));
}
uchar*
-RLCSA::readBWT()
+RLCSA::readBWT(pair_type range) const
{
- usint n = this->data_size + this->number_of_sequences;
+ if(isEmpty(range) || range.second >= this->data_size + this->number_of_sequences) { return 0; }
+
+ usint n = length(range);
uchar* bwt = new uchar[n];
memset(bwt, 0, n);
for(usint c = 0; c < CHARS; c++)
{
- RLEVector* vector = this->array[c];
- if(vector != 0)
+ if(this->array[c] != 0)
{
- bwt[vector->select(0)] = c;
- while(vector->hasNext()) { bwt[vector->selectNext()] = c; }
+ PsiVector::Iterator iter(*(this->array[c]));
+ usint pos = iter.valueAfter(range.first).first;
+ while(pos <= range.second)
+ {
+ bwt[pos - range.first] = c;
+ pos = iter.selectNext();
+ }
}
}
return bwt;
}
+usint
+RLCSA::countRuns() const
+{
+ usint runs = 0;
+ for(usint c = 0; c < CHARS; c++)
+ {
+ if(this->array[c] != 0)
+ {
+ PsiVector::Iterator iter(*(this->array[c]));
+ runs += iter.countRuns();
+ }
+ }
+
+ return runs;
+}
+
//--------------------------------------------------------------------------
-usint
-RLCSA::psi(usint index)
+PsiVector::Iterator**
+RLCSA::getIterators() const
{
- if(index >= this->data_size)
+ PsiVector::Iterator** iters = new PsiVector::Iterator*[CHARS];
+ for(usint i = 0; i < CHARS; i++)
{
- return this->data_size + this->number_of_sequences;
+ if(this->array[i] == 0) { iters[i] = 0; }
+ else { iters[i] = new PsiVector::Iterator(*(this->array[i])); }
}
+ return iters;
+}
- usint c = this->getCharacter(index);
- return this->array[c]->select(index - this->index_ranges[c].first);
+void
+RLCSA::deleteIterators(PsiVector::Iterator** iters) const
+{
+ if(iters == 0) { return; }
+ for(usint i = 0; i < CHARS; i++) { delete iters[i]; }
+ delete[] iters;
}
-pair_type
-RLCSA::psi(usint index, usint max_length)
+//--------------------------------------------------------------------------
+
+usint
+RLCSA::reportSize(bool print) const
{
- if(index >= this->data_size)
+ usint bytes = 0, temp = 0, bwt = 0;
+
+ for(usint c = 0; c < CHARS; c++)
+ {
+ if(this->array[c])
+ {
+ bytes += this->array[c]->reportSize();
+ bwt += this->array[c]->getCompressedSize();
+ }
+ }
+ bytes += sizeof(*this) + this->end_points->reportSize();
+ if(print)
{
- return pair_type(this->data_size + this->number_of_sequences, 0);
+ std::cout << "RLCSA: " << (bytes / (double)MEGABYTE) << " MB" << std::endl;
+ std::cout << " BWT only: " << (bwt / (double)MEGABYTE) << " MB" << std::endl;
}
- usint c = this->getCharacter(index);
- return this->array[c]->selectRun(index - this->index_ranges[c].first, max_length);
+ if(this->support_locate || this->support_display)
+ {
+ temp = this->sa_samples->reportSize();
+ if(print) { std::cout << "SA samples: " << (temp / (double)MEGABYTE) << " MB" << std::endl; }
+ bytes += temp;
+ }
+
+ if(print)
+ {
+ std::cout << "Total size: " << (bytes / (double)MEGABYTE) << " MB" << std::endl;
+ std::cout << std::endl;
+ }
+
+ return bytes;
+}
+
+void
+RLCSA::printInfo() const
+{
+ double megabytes = this->data_size / (double)MEGABYTE;
+
+ std::cout << "Sequences: " << this->number_of_sequences << std::endl;
+ std::cout << "Original size: " << megabytes << " MB" << std::endl;
+ std::cout << "Block size: " << (this->getBlockSize() * sizeof(usint)) << " bytes" << std::endl;
+ if(this->support_locate || this->support_display)
+ {
+ std::cout << "Sample rate: " << this->sample_rate << std::endl;
+ }
+ std::cout << std::endl;
}
//--------------------------------------------------------------------------
-pair_type
-RLCSA::getSequence(usint number)
+RLEVector*
+RLCSA::buildPLCP(usint block_size) const
{
- if(number >= this->number_of_sequences) { return EMPTY_PAIR; }
+ if(block_size < 2 * sizeof(usint) || block_size % sizeof(usint) != 0)
+ {
+ std::cerr << "PLCP: Block size must be a multiple of " << sizeof(usint) << " bytes!" << std::endl;
+ return 0;
+ }
- pair_type result;
- if(number == 0)
+ PsiVector::Iterator** iters = this->getIterators();
+
+ RLEEncoder plcp(block_size);
+ std::list<pair_type> matches;
+ pair_type prev_range = EMPTY_PAIR;
+ for(usint j = 0; j < this->number_of_sequences; j++)
{
- result.first = 0;
- result.second = this->end_points->select(number);
+ // Encode the padding as a single run.
+ pair_type seq_range = this->getSequenceRange(j);
+ if(j > 0 && prev_range.second + 1 < seq_range.first)
+ {
+ this->encodePLCPRun(plcp, prev_range.second + 1, seq_range.first, seq_range.first - prev_range.second - 2);
+ }
+ prev_range = seq_range;
+
+ usint maximal = seq_range.first;
+ usint x = this->sa_samples->inverseSA(seq_range.first), next_x;
+
+ // Invariant: x == inverseSA(i)
+ for(usint i = seq_range.first; i <= seq_range.second; i++, x = next_x)
+ {
+ usint c = this->getCharacter(x);
+
+ // T[i,n] is lexicographically the first suffix beginning with c.
+ if(x == this->index_ranges[c].first)
+ {
+ if(!matches.empty())
+ {
+ maximal = (*(matches.begin())).first;
+ matches.clear();
+ }
+ this->encodePLCPRun(plcp, maximal, i + 1, i - maximal);
+ maximal = i + 1;
+ next_x = this->psiUnsafe(x, c, *(iters[c])) - this->number_of_sequences;
+ continue;
+ }
+
+ // Process the previous left matches we are still following.
+ usint low = this->index_ranges[c].first + this->number_of_sequences;
+ usint start, stop; start = stop = maximal;
+ for(std::list<pair_type>::iterator iter = matches.begin(); iter != matches.end();)
+ {
+ pair_type match = *iter;
+ if(match.second < low) // These no longer match the current character.
+ {
+ if(match.first < start) { start = match.first; } // Form a single run from then.
+ iter = matches.erase(iter);
+ }
+ else
+ {
+ if(match.first < stop) { stop = match.first; } // End of the run to be encoded.
+ (*iter).second = this->psiUnsafe(match.second - this->number_of_sequences, c, *(iters[c]));
+ ++iter;
+ }
+ }
+ if(start < stop) { this->encodePLCPRun(plcp, start, stop, i - start); }
+
+ // If PLCP[i] is minimal, we add the left match to the list.
+ // We add pairs of type (j, inverseSA(j) + n_of_s) as inverseSA is < 0 for end markers.
+ usint next_y = this->psiUnsafe(x - 1, c, *(iters[c]));
+ next_x = this->psiUnsafeNext(c, *(iters[c])) - this->number_of_sequences;
+ if(next_y != next_x + this->number_of_sequences - 1)
+ {
+ matches.push_back(pair_type(maximal, next_y));
+ maximal = i + 1;
+ }
+ }
+
+ if(!matches.empty())
+ {
+ maximal = (*(matches.begin())).first;
+ matches.clear();
+ }
+ if(maximal <= seq_range.second)
+ {
+ this->encodePLCPRun(plcp, maximal, seq_range.second + 1, seq_range.second + 1 - maximal);
+ }
}
- else
+
+ this->deleteIterators(iters);
+
+ return new RLEVector(plcp, 2 * (this->end_points->getSize() + 1));
+}
+
+usint
+RLCSA::sampleLCP(usint sample_rate, pair_type*& sampled_values, bool report) const
+{
+ if(sample_rate == 0)
{
- if(this->sa_samples == 0) { return EMPTY_PAIR; }
- usint d = this->sa_samples->getSampleRate();
- result.first = d * ((this->end_points->select(number - 1) / d) + 1);
- result.second = this->end_points->selectNext();
+ sample_rate = this->data_size + 1;
}
- return result;
+ PsiVector::Iterator** iters = this->getIterators();
+
+ usint runs = this->countRuns();
+ usint samples = 0, minimal_sum = 0, nonstrict_minimal_sum = 0, minimal_samples = 0;
+ usint max_samples = runs + (this->data_size - runs) / sample_rate;
+ sampled_values = new pair_type[max_samples];
+
+ std::list<Triple> matches;
+ for(usint j = 0; j < this->number_of_sequences; j++)
+ {
+ pair_type seq_range = this->getSequenceRange(j);
+ usint first_sample = samples; // First minimal sample of the current sequence.
+ usint minimal = seq_range.second + 1; // Last minimal position.
+ usint i, x = this->sa_samples->inverseSA(seq_range.first), next_x, prev_x = 0;
+
+ // Invariant: x == inverseSA(i)
+ for(i = seq_range.first; i <= seq_range.second; i++, x = next_x)
+ {
+ usint c = this->getCharacter(x);
+
+ // T[i,n] is lexicographically the first suffix beginning with c.
+ if(x == this->index_ranges[c].first)
+ {
+ if(!matches.empty())
+ {
+ for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end(); ++iter)
+ {
+ nonstrict_minimal_sum += i - (*iter).first;
+ }
+ matches.clear();
+ }
+
+ // Minimal sample: SA[x] = i, LCP[x] = 0
+ sampled_values[samples] = pair_type(x, 0); minimal = i;
+ samples++; minimal_samples++;
+ next_x = this->psiUnsafe(x, c, *(iters[c])) - this->number_of_sequences; prev_x = x;
+ continue;
+ }
+
+ // Process the previous left matches we are still following.
+ usint low = this->index_ranges[c].first + this->number_of_sequences;
+ pair_type potential_sample(this->data_size, this->data_size);
+ for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end();)
+ {
+ Triple match = *iter;
+ if(match.second < low) // These no longer match the current character.
+ {
+ if(potential_sample.first != this->data_size) { nonstrict_minimal_sum += potential_sample.second; }
+
+ // Potential minimal sample: SA[match.third] = match.first, LCP[match.third] = i - match.first
+ potential_sample = pair_type(match.third, i - match.first);
+ iter = matches.erase(iter);
+ }
+ else
+ {
+ (*iter).second = this->psiUnsafe(match.second - this->number_of_sequences, c, *(iters[c]));
+ ++iter;
+ }
+ }
+ if(potential_sample.first != this->data_size)
+ {
+ // Last potential sample is minimal.
+ sampled_values[samples] = potential_sample;
+ samples++; minimal_sum += potential_sample.second; minimal_samples++;
+ }
+
+ // If PLCP[i] is minimal, we add the left match to the list.
+ // We add triples of type (j, inverseSA(j) + n_of_s, x) as inverseSA is < 0 for end markers.
+ // x is the original minimal position.
+ usint next_y = this->psiUnsafe(x - 1, c, *(iters[c]));
+ next_x = this->psiUnsafeNext(c, *(iters[c]));
+ if(next_y != next_x - 1 || next_x < this->number_of_sequences)
+ {
+ matches.push_back(Triple(i, next_y, x)); minimal = i;
+ }
+ next_x -= this->number_of_sequences; prev_x = x;
+ }
+
+ // Are we still following something?
+ if(!matches.empty())
+ {
+ pair_type potential_sample(this->data_size, this->data_size);
+ for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end(); ++iter)
+ {
+ Triple match = *iter;
+ if(potential_sample.first != this->data_size) { nonstrict_minimal_sum += potential_sample.second; }
+
+ // Potential minimal sample: SA[match.third] = match.first, LCP[match.third] = i - match.first
+ potential_sample = pair_type(match.third, i - match.first);
+ }
+ matches.clear();
+ if(potential_sample.first != this->data_size)
+ {
+ // Last potential sample is minimal.
+ sampled_values[samples] = potential_sample;
+ samples++; minimal_sum += potential_sample.second; minimal_samples++;
+ }
+ }
+
+ // Add the non-minimal samples.
+ if(sample_rate <= this->data_size)
+ {
+ usint last_sample = samples - 1;
+ i = seq_range.first; x = this->sa_samples->inverseSA(seq_range.first);
+ for(usint current_sample = first_sample; current_sample <= last_sample; current_sample++)
+ {
+ // Find the next minimal sample and add nonminimal samples if needed.
+ pair_type first_nonminimal(i + sample_rate - 1, samples);
+ pair_type next_nonminimal = first_nonminimal;
+ while(x != sampled_values[current_sample].first)
+ {
+ if(i == next_nonminimal.first)
+ {
+ sampled_values[samples] = pair_type(x, 0); samples++;
+ next_nonminimal.first += sample_rate; next_nonminimal.second++;
+ }
+ i++; x = this->psi(x, iters) - this->number_of_sequences;
+ }
+
+ // Reduce the nonminimal samples to the current minimal sample.
+ for(next_nonminimal = first_nonminimal; next_nonminimal.second < samples; next_nonminimal.first += sample_rate, next_nonminimal.second++)
+ {
+ sampled_values[next_nonminimal.second].second = sampled_values[current_sample].second + i - next_nonminimal.first;
+ }
+ i++; x = this->psi(x, iters) - this->number_of_sequences;
+ }
+ }
+ }
+
+ std::sort(sampled_values, sampled_values + samples);
+
+ if(report)
+ {
+ std::cout << "Samples: " << samples << " (total) / " << minimal_samples << " (minimal)" << std::endl;
+ std::cout << "Upper bounds: " << max_samples << " (total) / " << runs << " (minimal)" << std::endl;
+ std::cout << "Sum of minimal samples: " << (minimal_sum + nonstrict_minimal_sum) << " (total) / " << minimal_sum << " (strict)" << std::endl;
+ std::cout << std::endl;
+ }
+
+ this->deleteIterators(iters);
+
+ return samples;
}
usint
-RLCSA::reportSize(bool print)
+RLCSA::lcp(usint index, const LCPSamples& lcp_samples) const
{
- usint bytes = 0, temp = 0;
+ if(index >= this->data_size) { return this->data_size; }
- for(usint c = 0; c < CHARS; c++)
+ usint offset = 0;
+ index += this->number_of_sequences;
+ while(true)
{
- if(this->array[c]) { temp += this->array[c]->reportSize(); }
+ pair_type next_sample = lcp_samples.getFirstSampleAfter(index - this->number_of_sequences);
+ next_sample.first += this->number_of_sequences;
+ if(next_sample.first == index)
+ {
+ return lcp_samples.getSample(next_sample.second) + offset;
+ }
+ index = this->psi(index - this->number_of_sequences);
+ offset++;
}
- if(print) { std::cout << "Psi array: " << temp << std::endl; }
- bytes += temp;
+}
- if(this->support_locate || this->support_display)
+usint
+RLCSA::lcp(usint index, const LCPSamples& lcp_samples, const RLEVector& plcp) const
+{
+ if(index >= this->data_size) { return this->data_size; }
+
+ usint offset = 0;
+ while(true)
{
- temp = this->sa_samples->reportSize();
- if(print) { std::cout << "SA samples: " << temp << std::endl; }
- bytes += temp;
+ usint c = this->getCharacter(index);
+ PsiVector::Iterator iter(*(this->array[c]));
+
+ if(index == this->index_ranges[c].first)
+ {
+ return lcp_samples.getSampleAtPosition(index) + offset;
+ }
+ usint next_y = this->psiUnsafe(index - 1, c, iter);
+ usint next_x = this->psiUnsafeNext(c, iter);
+ if(next_y != next_x - 1 || next_x < this->number_of_sequences)
+ {
+ pair_type sample = lcp_samples.getFirstSampleAfter(index);
+ if(sample.first == index)
+ {
+ return lcp_samples.getSample(sample.second) + offset;
+ }
+ }
+
+ pair_type next_sample = this->sa_samples->getFirstSampleAfter(index);
+ if(next_sample.first == index)
+ {
+ index = this->sa_samples->getSample(next_sample.second) - offset;
+ RLEVector::Iterator plcp_iter(plcp);
+ return plcp_iter.select(index) - 2 * index;
+ }
+ index = next_x - this->number_of_sequences;
+ offset++;
+ }
+}
+
+usint
+RLCSA::lcpDirect(usint index) const
+{
+ if(index >= this->data_size || index == 0)
+ {
+ return 0;
}
- temp = sizeof(*this) + this->end_points->reportSize();
- if(print) { std::cout << "RLCSA overhead: " << temp << std::endl; }
- bytes += temp;
+ usint match = index - 1;
+ usint value = 0;
- if(print)
+ usint c = this->getCharacter(index);
+ usint low = this->index_ranges[c].first;
+ while(match >= low)
{
- std::cout << "Total size: " << bytes << std::endl;
- std::cout << std::endl;
+ PsiVector::Iterator iter(*(this->array[c]));
+ match = this->psiUnsafe(match, c, iter);
+ index = this->psiUnsafe(index, c, iter);
+ value++;
+
+ if(match < this->number_of_sequences || index < this->number_of_sequences)
+ {
+ break;
+ }
+ match -= this->number_of_sequences;
+ index -= this->number_of_sequences;
+
+ c = this->getCharacter(index);
+ low = this->index_ranges[c].first;
}
- return bytes;
+ return value;
+}
+
+//--------------------------------------------------------------------------
+
+void
+RLCSA::mergeEndPoints(RLCSA& index, RLCSA& increment)
+{
+ DeltaEncoder* endings = new DeltaEncoder(RLCSA::ENDPOINT_BLOCK_SIZE);
+
+ DeltaVector::Iterator index_iter(*(index.end_points));
+ DeltaVector::Iterator increment_iter(*(increment.end_points));
+
+ endings->setBit(index_iter.select(0));
+ for(usint i = 1; i < index.number_of_sequences; i++)
+ {
+ endings->setBit(index_iter.selectNext());
+ }
+ usint sum = index.end_points->getSize();
+ delete index.end_points; index.end_points = 0;
+
+ endings->setBit(sum + increment_iter.select(0));
+ for(usint i = 1; i < increment.number_of_sequences; i++)
+ {
+ endings->setBit(sum + increment_iter.selectNext());
+ }
+ sum += increment.end_points->getSize();
+ delete increment.end_points; increment.end_points = 0;
+
+ this->end_points = new DeltaVector(*endings, sum);
+ delete endings;
+}
+
+
+void
+RLCSA::mergeSamples(RLCSA& index, RLCSA& increment, usint* positions)
+{
+ if(this->support_locate || this->support_display)
+ {
+ positions += increment.number_of_sequences;
+ this->sa_samples = new SASamples(*(index.sa_samples), *(increment.sa_samples), positions, increment.data_size, this->number_of_sequences);
+ }
}
//--------------------------------------------------------------------------
#define RLCSA_H
#include <fstream>
-#include <cstring> // defines std::memset, added by Kim
+#include <vector>
+
+#include "bits/deltavector.h"
+#include "bits/rlevector.h"
+#include "bits/nibblevector.h"
-#include "bits/vectors.h"
#include "sasamples.h"
+#include "lcpsamples.h"
#include "misc/parameters.h"
const std::string ARRAY_EXTENSION = ".rlcsa.array";
const std::string SA_SAMPLES_EXTENSION = ".rlcsa.sa_samples";
const std::string PARAMETERS_EXTENSION = ".rlcsa.parameters";
+const std::string LCP_SAMPLES_EXTENSION = ".lcp_samples";
+const std::string PLCP_EXTENSION = ".plcp";
const parameter_type RLCSA_BLOCK_SIZE = parameter_type("RLCSA_BLOCK_SIZE", 32);
const parameter_type SAMPLE_RATE = parameter_type("SAMPLE_RATE", 512);
-const parameter_type SUPPORT_LOCATE = parameter_type("SUPPORT_LOCATE", 0);
-const parameter_type SUPPORT_DISPLAY = parameter_type("SUPPORT_DISPLAY", 0);
+const parameter_type SUPPORT_LOCATE = parameter_type("SUPPORT_LOCATE", 1);
+const parameter_type SUPPORT_DISPLAY = parameter_type("SUPPORT_DISPLAY", 1);
-struct LocateItem
-{
- usint value;
- usint offset;
- bool found;
-};
-
+#ifdef USE_NIBBLE_VECTORS
+typedef NibbleVector PsiVector;
+#else
+typedef RLEVector PsiVector;
+#endif
class RLCSA
{
public:
+
+//--------------------------------------------------------------------------
+// CONSTRUCTION
+//--------------------------------------------------------------------------
+
const static usint ENDPOINT_BLOCK_SIZE = 16;
RLCSA(const std::string& base_name, bool print = false);
/*
If multiple_sequences is true, each 0 is treated as a end of sequence marker.
There must be nonzero characters between the 0s. data[bytes - 1] must also be 0.
+ FIXME Crashes if bytes >= 2 GB.
*/
RLCSA(uchar* data, usint bytes, usint block_size, usint sa_sample_rate, bool multiple_sequences, bool delete_data);
// Destroys contents of index and increment.
- RLCSA(RLCSA& index, RLCSA& increment, usint* positions, usint block_size);
+ // threads has no effect unless MULTITHREAD_SUPPORT is defined
+ RLCSA(RLCSA& index, RLCSA& increment, usint* positions, usint block_size, usint threads = 1);
~RLCSA();
- void writeTo(const std::string& base_name);
+ void writeTo(const std::string& base_name) const;
+
+ bool isOk() const;
- bool isOk();
+//--------------------------------------------------------------------------
+// QUERIES
+//--------------------------------------------------------------------------
// Returns the closed interval of suffix array containing the matches.
- pair_type count(const std::string& pattern);
+ pair_type count(const std::string& pattern) const;
- void reportPositions(uchar* data, usint length, usint* positions);
+ // Used when merging CSAs.
+ void reportPositions(uchar* data, usint length, usint* positions) const;
// Returns SA[range]. User must free the buffer. Latter version uses buffer provided by the user.
- LocateItem* locate(pair_type range);
- LocateItem* locate(pair_type range, LocateItem* data);
+ usint* locate(pair_type range) const;
+ usint* locate(pair_type range, usint* data) const;
// Returns SA[index].
- usint locate(usint index);
+ usint locate(usint index) const;
+
+ // Returns T^{sequence}[range]. User must free the buffer.
+ // Third version uses buffer provided by the user.
+ uchar* display(usint sequence) const;
+ uchar* display(usint sequence, pair_type range) const;
+ uchar* display(usint sequence, pair_type range, uchar* data) const;
+
+ // Displays the intersection of T[position - context, position + len + context - 1]
+ // and T^{getSequenceForPosition(position)}.
+ // This is intended for displaying an occurrence of a pattern of length 'len'
+ // with 'context' extra characters on both sides.
+ // The actual length of the returned string is written into result_length.
+ uchar* display(usint position, usint len, usint context, usint& result_length) const;
- // Returns Ti[range]. User must free the buffer. Latter version uses buffer provided by the user.
- uchar* display(usint sequence, pair_type range);
- uchar* display(usint sequence, pair_type range, uchar* data);
+ // Returns the actual length of the prefix. User must provide the buffer.
+ usint displayPrefix(usint sequence, usint len, uchar* data) const;
- // Writes the Psi array into given file. End of sequence markers are not written.
- void decompressInto(std::ofstream& psi_file);
+ // Returns at most max_len characters starting from T[SA[index]].
+ // User must provide the buffer. Returns the number of characters in buffer.
+ usint displayFromPosition(usint index, usint max_len, uchar* data) const;
+
+ // Get the range of SA values for the sequence identified by
+ // a sequence number or a SA value.
+ pair_type getSequenceRange(usint number) const;
+ pair_type getSequenceRangeForPosition(usint value) const;
+
+ // Get the sequence number for given SA value(s).
+ // The returned array is the same as the parameter.
+ usint getSequenceForPosition(usint value) const;
+ usint* getSequenceForPosition(usint* value, usint length) const;
+
+ // Changes SA value to (sequence, offset).
+ pair_type getRelativePosition(usint value) const;
// Returns the BWT of the collection including end of sequence markers.
- uchar* readBWT();
+ uchar* readBWT() const;
+ uchar* readBWT(pair_type range) const;
+
+ // Returns the number of equal letter runs in the BWT. Runs consisting of end markers are ignored.
+ usint countRuns() const;
- // Return value includes the implicit end of sequence markers. To get suffix array indexes,
- // subtract getNumberOfSequences() from the value.
- usint psi(usint index);
- pair_type psi(usint index, usint max_length); // This version returns a run.
+//--------------------------------------------------------------------------
+// BASIC OPERATIONS
+//--------------------------------------------------------------------------
- inline bool supportsLocate() { return this->support_locate; }
- inline bool supportsDisplay() { return this->support_display; }
- inline usint getSize() { return this->data_size; }
- inline usint getNumberOfSequences() { return this->number_of_sequences; }
+ // The return values of these functions include the implicit end markers.
+ // To get SA indexes, subtract getNumberOfSequences() from the value.
- pair_type getSequence(usint number);
+ inline usint psi(usint index) const
+ {
+ if(index >= this->data_size)
+ {
+ return this->data_size + this->number_of_sequences;
+ }
+
+ usint c = this->getCharacter(index);
+ return this->psiUnsafe(index, c);
+ }
+
+ // This version returns a run.
+ inline pair_type psi(usint index, usint max_length) const
+ {
+ if(index >= this->data_size)
+ {
+ return pair_type(this->data_size + this->number_of_sequences, 0);
+ }
+
+ usint c = this->getCharacter(index);
+ PsiVector::Iterator iter(*(this->array[c]));
+ return iter.selectRun(index - this->index_ranges[c].first, max_length);
+ }
+
+ inline usint C(usint c) const
+ {
+ if(c >= CHARS)
+ return this->data_size + this->number_of_sequences;
+ return this->index_ranges[c].first + this->number_of_sequences - 1;
+ }
+
+ inline usint LF(usint index, usint c) const
+ {
+ if(c >= CHARS)
+ {
+ return this->data_size + this->number_of_sequences;
+ }
+ if(this->array[c] == 0)
+ {
+ if(c < this->text_chars[0]) { return this->number_of_sequences - 1; }
+ return this->index_ranges[c].first + this->number_of_sequences - 1;
+ }
+ index += this->number_of_sequences;
+
+ PsiVector::Iterator iter(*(this->array[c]));
+ return this->LF(index, c, iter);
+ }
+
+//--------------------------------------------------------------------------
+// REPORTING
+//--------------------------------------------------------------------------
+
+ inline bool supportsLocate() const { return this->support_locate; }
+ inline bool supportsDisplay() const { return this->support_display; }
+ inline usint getSize() const { return this->data_size; }
+ inline usint getNumberOfSequences() const { return this->number_of_sequences; }
+ inline usint getBlockSize() const { return this->array[this->text_chars[0]]->getBlockSize(); }
// Returns the size of the data structure.
- usint reportSize(bool print = false);
+ usint reportSize(bool print = false) const;
+
+ void printInfo() const;
+
+//--------------------------------------------------------------------------
+// LCP EXPERIMENTS
+//--------------------------------------------------------------------------
+
+ // Optimized version:
+ // - Interleaves main loop with computing irreducible values.
+ // - Encodes maximal runs from a true local maximum to a true local minimum.
+ RLEVector* buildPLCP(usint block_size) const;
+
+ // Returns the number of samples. sampled_values will be a pointer to the samples.
+ usint sampleLCP(usint sample_rate, pair_type*& sampled_values, bool report = false) const;
+
+ usint lcp(usint index, const LCPSamples& lcp_samples) const;
+ usint lcp(usint index, const LCPSamples& lcp_samples, const RLEVector& plcp) const;
+
+ // Calculate LCP[index] directly.
+ usint lcpDirect(usint index) const;
+
+ // Writes PLCP[start] to PLCP[stop - 1].
+ inline void encodePLCPRun(RLEEncoder& plcp, usint start, usint stop, usint first_val) const
+ {
+ plcp.setRun(2 * start + first_val, stop - start);
+// std::cerr << "(" << start << ", " << stop << ", " << first_val << ")" << std::endl;
+ }
+
+//--------------------------------------------------------------------------
+// INTERNAL VARIABLES
+//--------------------------------------------------------------------------
private:
bool ok;
- RLEVector* array[CHARS];
+ PsiVector* array[CHARS];
SASamples* sa_samples;
pair_type index_ranges[CHARS];
usint number_of_sequences;
DeltaVector* end_points;
- void locateUnsafe(pair_type range, LocateItem* data);
- bool processRun(pair_type run, LocateItem* data);
- void displayUnsafe(pair_type range, uchar* data);
+//--------------------------------------------------------------------------
+// INTERNAL VERSIONS OF QUERIES
+//--------------------------------------------------------------------------
+
+ void locateUnsafe(pair_type range, usint* data) const;
+ bool processRun(pair_type run, usint* data, usint* offsets, bool* finished, PsiVector::Iterator** iters) const;
+ void displayUnsafe(pair_type range, uchar* data) const;
+
+//--------------------------------------------------------------------------
+// INTERNAL VERSIONS OF BASIC OPERATIONS
+//--------------------------------------------------------------------------
+
+ inline usint psi(usint index, PsiVector::Iterator** iters) const
+ {
+ usint c = this->getCharacter(index);
+ return iters[c]->select(index - this->index_ranges[c].first);
+ }
- inline usint getCharacter(usint pos)
+ inline pair_type psi(usint index, usint max_length, PsiVector::Iterator** iters) const
{
- usint* curr = &(this->text_chars[this->index_pointers[pos / this->index_rate]]);
+ usint c = this->getCharacter(index);
+ return iters[c]->selectRun(index - this->index_ranges[c].first, max_length);
+ }
+
+ // Returns psi(index), assuming the suffix of rank index begins with c.
+ inline usint psiUnsafe(usint index, usint c) const
+ {
+ PsiVector::Iterator iter(*(this->array[c]));
+ return this->psiUnsafe(index, c, iter);
+ }
+
+ // As above, but with a given iterator.
+ inline usint psiUnsafe(usint index, usint c, PsiVector::Iterator& iter) const
+ {
+ return iter.select(index - this->index_ranges[c].first);
+ }
+
+ // As above, but returns the next value of psi.
+ inline usint psiUnsafeNext(usint c, PsiVector::Iterator& iter) const
+ {
+ return iter.selectNext();
+ }
+
+ inline pair_type backwardSearchStep(pair_type range, usint c) const
+ {
+ if(this->array[c] == 0) { return EMPTY_PAIR; }
+ PsiVector::Iterator iter(*(this->array[c]));
+
+ usint start = this->index_ranges[c].first + this->number_of_sequences - 1;
+ range.first = start + iter.rank(range.first, true);
+ range.second = start + iter.rank(range.second);
+
+ return range;
+ }
+
+ inline usint LF(usint index, usint c, PsiVector::Iterator& iter) const
+ {
+ return this->index_ranges[c].first + this->number_of_sequences + iter.rank(index) - 1;
+ }
+
+//--------------------------------------------------------------------------
+// INTERNAL STUFF
+//--------------------------------------------------------------------------
+
+ // Creates an array of iterators for every vector in this->array.
+ PsiVector::Iterator** getIterators() const;
+ void deleteIterators(PsiVector::Iterator** iters) const;
+
+ inline usint getCharacter(usint pos) const
+ {
+ const usint* curr = &(this->text_chars[this->index_pointers[pos / this->index_rate]]);
while(pos > this->index_ranges[*curr].second) { curr++; }
return *curr;
}
+ void mergeEndPoints(RLCSA& index, RLCSA& increment);
+ void mergeSamples(RLCSA& index, RLCSA& increment, usint* positions);
+
void buildCharIndexes(usint* distribution);
+
+ // These are not allowed.
+ RLCSA();
+ RLCSA(const RLCSA&);
+ RLCSA& operator = (const RLCSA&);
};
+#include <algorithm>
+#include <cstdlib>
+#include <cstring>
#include <iostream>
#include "rlcsa_builder.h"
+#include "misc/utils.h"
+
+#ifdef MULTITHREAD_SUPPORT
+#include <omp.h>
+#endif
namespace CSA
{
-RLCSABuilder::RLCSABuilder(usint _block_size, usint _sample_rate, usint _buffer_size) :
+RLCSABuilder::RLCSABuilder(usint _block_size, usint _sample_rate, usint _buffer_size, usint _threads) :
block_size(_block_size), sample_rate(_sample_rate), buffer_size(_buffer_size),
+ threads(_threads),
buffer(0)
{
this->reset();
+ this->build_time = this->search_time = this->sort_time = this->merge_time = 0.0;
}
RLCSABuilder::~RLCSABuilder()
if(this->buffer == 0)
{
- clock_t start = clock();
+ double start = readTimer();
RLCSA* temp = new RLCSA((uchar*)sequence, length, this->block_size, this->sample_rate, false, false);
- this->build_time += clock() - start;
+ this->build_time += readTimer() - start;
this->addRLCSA(temp, (uchar*)sequence, length + 1, delete_sequence);
return;
}
}
else
{
- this->flush();
- this->buffer = new uchar[this->buffer_size];
+ if(this->chars > 0)
+ {
+ this->flush();
+ this->buffer = new uchar[this->buffer_size];
+ }
if(length >= this->buffer_size - 1)
{
- clock_t start = clock();
+ double start = readTimer();
RLCSA* temp = new RLCSA((uchar*)sequence, length, this->block_size, this->sample_rate, false, false);
- this->build_time += clock() - start;
+ this->build_time += readTimer() - start;
this->addRLCSA(temp, (uchar*)sequence, length + 1, delete_sequence);
}
else
}
}
+void
+RLCSABuilder::insertFromFile(const std::string& base_name)
+{
+ if(!this->ok) { return; }
+
+ if(this->buffer != 0 && this->chars > 0)
+ {
+ this->flush();
+ this->buffer = new uchar[this->buffer_size];
+ }
+
+ std::ifstream input(base_name.c_str(), std::ios_base::binary);
+ if(!input) { return; }
+ RLCSA* increment = new RLCSA(base_name);
+ usint data_size = increment->getSize() + increment->getNumberOfSequences();
+ uchar* data = new uchar[data_size];
+ input.read((char*)data, data_size);
+ input.close();
+
+ this->addRLCSA(increment, data, data_size, true);
+}
+
RLCSA*
RLCSABuilder::getRLCSA()
{
double
RLCSABuilder::getBuildTime()
{
- return this->build_time / (double)CLOCKS_PER_SEC;
+ return this->build_time;
}
double
RLCSABuilder::getSearchTime()
{
- return this->search_time / (double)CLOCKS_PER_SEC;
+ return this->search_time;
+}
+
+double
+RLCSABuilder::getSortTime()
+{
+ return this->sort_time;
}
double
RLCSABuilder::getMergeTime()
{
- return this->merge_time / (double)CLOCKS_PER_SEC;
+ return this->merge_time;
}
//--------------------------------------------------------------------------
void
RLCSABuilder::flush()
{
- clock_t start = clock();
+ double start = readTimer();
RLCSA* temp = new RLCSA(this->buffer, this->chars, this->block_size, this->sample_rate, true, (this->index == 0));
- this->build_time += clock() - start;
- this->addRLCSA(temp, this->buffer, this->chars, true);
+ this->build_time += readTimer() - start;
+ this->addRLCSA(temp, this->buffer, this->chars, (this->index != 0));
this->buffer = 0; this->chars = 0;
}
{
if(this->index != 0)
{
- clock_t start = clock();
+ double start = readTimer();
- usint* positions = new usint[length];
- usint begin = 0;
+ usint sequences = increment->getNumberOfSequences();
+ usint* end_markers = new usint[sequences];
+ usint curr = 0;
for(usint i = 0; i < length - 1; i++)
{
- if(sequence[i] == 0)
+ if(sequence[i] == 0) { end_markers[curr++] = i; }
+ }
+ end_markers[sequences - 1] = length - 1;
+
+ usint* positions = new usint[length]; usint begin;
+ #ifdef MULTITHREAD_SUPPORT
+ usint chunk = std::max((usint)1, sequences / (8 * this->threads));
+ omp_set_num_threads(this->threads);
+ #pragma omp parallel private(begin)
+ {
+ #pragma omp for schedule(dynamic, chunk)
+ for(sint i = 0; i < (sint)sequences; i++)
{
- this->index->reportPositions(&(sequence[begin]), i - begin, &(positions[begin]));
- begin = i + 1;
+ if(i > 0) { begin = end_markers[i - 1] + 1; } else { begin = 0; }
+ this->index->reportPositions(sequence + begin, end_markers[i] - begin, positions + begin);
}
}
- this->index->reportPositions(&(sequence[begin]), length - 1 - begin, &(positions[begin]));
+ #else
+ for(sint i = 0; i < (sint)sequences; i++)
+ {
+ if(i > 0) { begin = end_markers[i - 1] + 1; } else { begin = 0; }
+ this->index->reportPositions(sequence + begin, end_markers[i] - begin, positions + begin);
+ }
+ #endif
+ delete[] end_markers;
+ if(delete_sequence) { delete[] sequence; }
+
+ double mark = readTimer();
+ this->search_time += mark - start;
+ #ifdef MULTITHREAD_SUPPORT
+ omp_set_num_threads(this->threads);
+ #endif
std::sort(positions, positions + length);
for(usint i = 0; i < length; i++)
{
positions[i] += i + 1; // +1 because the insertion will be after positions[i]
}
- if(delete_sequence) { delete[] sequence; }
- clock_t mark = clock();
- this->search_time += mark - start;
+ double sort = readTimer();
+ this->sort_time += sort - mark;
- RLCSA* merged = new RLCSA(*(this->index), *increment, positions, this->block_size);
+ RLCSA* merged = new RLCSA(*(this->index), *increment, positions, this->block_size, this->threads);
delete[] positions;
delete this->index;
delete increment;
this->index = merged;
- this->merge_time += clock() - mark;
+ this->merge_time += readTimer() - sort;
}
else
{
+ if(delete_sequence) { delete[] sequence; }
this->index = increment;
}
#ifndef RLCSA_BUILDER_H
#define RLCSA_BUILDER_H
-#include <cstdlib>
#include "rlcsa.h"
class RLCSABuilder
{
public:
- RLCSABuilder(usint _block_size, usint _sample_rate, usint _buffer_size);
+ RLCSABuilder(usint _block_size, usint _sample_rate, usint _buffer_size, usint _threads = 1);
~RLCSABuilder();
void insertSequence(char* sequence, usint length, bool delete_sequence);
+ // Use this if you have already built an index for the file.
+ void insertFromFile(const std::string& base_name);
+
// User must free the index. Builder no longer contains it.
RLCSA* getRLCSA();
// These times are not reset with the rest of the builder.
double getBuildTime();
double getSearchTime();
+ double getSortTime();
double getMergeTime();
private:
usint sample_rate;
usint buffer_size;
+ usint threads;
+
uchar* buffer;
usint chars;
bool ok;
- clock_t build_time;
- clock_t search_time;
- clock_t merge_time;
+ double build_time;
+ double search_time;
+ double sort_time;
+ double merge_time;
void flush();
void addRLCSA(RLCSA* increment, uchar* sequence, usint length, bool delete_sequence);
void reset();
+
+ // These are not allowed.
+ RLCSABuilder();
+ RLCSABuilder(const RLCSABuilder&);
+ RLCSABuilder& operator = (const RLCSABuilder&);
};
--- /dev/null
+#include <iostream>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+using namespace CSA;
+
+
+enum mode_type { COUNT, DISPLAY, CONTEXT };
+
+
+void printUsage()
+{
+ std::cout << "Usage: rlcsa_grep [-c] pattern base_name" << std::endl;
+ std::cout << " -c print the number of matching lines" << std::endl;
+ std::cout << " -NUM display NUM characters of leading and trailing context instead of" << std::endl;
+ std::cout << " the entire line" << std::endl;
+}
+
+
+int main(int argc, char** argv)
+{
+ int base_arg = 2, pattern_arg = 1;
+ mode_type mode = DISPLAY;
+ usint context = 0;
+
+ if(argc < base_arg + 1)
+ {
+ printUsage();
+ return 1;
+ }
+
+ if(argv[1][0] == '-')
+ {
+ base_arg++; pattern_arg++;
+ if(std::string("-c").compare(argv[1]) == 0)
+ {
+ mode = COUNT;
+ }
+ else
+ {
+ mode = CONTEXT;
+ context = atoi(&(argv[1][1]));
+ }
+ if(argc < base_arg + 1)
+ {
+ printUsage();
+ return 2;
+ }
+ }
+
+ RLCSA rlcsa(argv[base_arg], false);
+ if(!rlcsa.isOk())
+ {
+ return 3;
+ }
+
+ usint len = std::string(argv[pattern_arg]).length();
+ pair_type result_range = rlcsa.count(argv[pattern_arg]);
+ usint occurrences = length(result_range);
+ if(occurrences == 0)
+ {
+ if(mode == COUNT)
+ {
+ std::cout << 0 << std::endl;
+ }
+ return 0;
+ }
+
+ usint last_row = 0;
+ usint* results = rlcsa.locate(result_range);
+ if(mode != CONTEXT)
+ {
+ rlcsa.getSequenceForPosition(results, occurrences);
+ }
+ std::sort(results, results + occurrences);
+ if(mode != CONTEXT)
+ {
+ for(usint i = 1; i < occurrences; i++)
+ {
+ if(results[i] != results[last_row])
+ {
+ last_row++; results[last_row] = results[i];
+ }
+ }
+ }
+
+ if(mode == COUNT)
+ {
+ std::cout << (last_row + 1) << std::endl;
+ }
+ else if(mode == DISPLAY)
+ {
+ for(usint i = 0; i <= last_row; i++)
+ {
+ uchar* row = rlcsa.display(results[i]);
+ std::cout.write((char*)row, length(rlcsa.getSequenceRange(results[i])));
+ std::cout << std::endl;
+ delete[] row;
+ }
+ }
+ else if(mode == CONTEXT)
+ {
+ usint result_length = 0;
+ for(usint i = 0; i < occurrences; i++)
+ {
+ uchar* text = rlcsa.display(results[i], len, context, result_length);
+ std::cout.write((char*)text, result_length);
+ std::cout << std::endl;
+ delete[] text;
+ }
+ }
+
+ delete[] results;
+ return 0;
+}
--- /dev/null
+#include <cstdlib>
+#include <ctime>
+#include <iostream>
+#include <vector>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+#ifdef MULTITHREAD_SUPPORT
+#include <omp.h>
+#endif
+
+
+using namespace CSA;
+
+
+const int MAX_THREADS = 64;
+const usint MAX_OCCURRENCES = 100000;
+
+
+int main(int argc, char** argv)
+{
+ std::cout << "RLCSA test" << std::endl;
+ if(argc < 2)
+ {
+ std::cout << "Usage: rlcsa_test basename [patterns [threads]]" << std::endl;
+ return 1;
+ }
+
+ std::cout << "Base name: " << argv[1] << std::endl;
+ if(argc > 2) { std::cout << "Patterns: " << argv[2] << std::endl; }
+ sint threads = 1;
+ if(argc > 3)
+ {
+ threads = std::min(MAX_THREADS, std::max(atoi(argv[3]), 1));
+ }
+ std::cout << "Threads: " << threads << std::endl;
+ std::cout << std::endl;
+
+ const RLCSA rlcsa(argv[1], false);
+ rlcsa.printInfo();
+ rlcsa.reportSize(true);
+
+ if(argc < 3) { return 0; }
+ std::ifstream patterns(argv[2], std::ios_base::binary);
+ if(!patterns)
+ {
+ std::cerr << "Error opening pattern file!" << std::endl;
+ return 2;
+ }
+
+ std::vector<std::string> rows;
+ readRows(patterns, rows, true);
+
+ usint total = 0, total_size = 0, ignored = 0, found = 0;
+ sint i, n = rows.size();
+ pair_type result;
+ usint* matches;
+
+ double start = readTimer();
+ #ifdef MULTITHREAD_SUPPORT
+ usint occurrences;
+ omp_set_num_threads(threads);
+ #pragma omp parallel private(result, matches, occurrences)
+ {
+ #pragma omp for schedule(dynamic, 1)
+ for(i = 0; i < n; i++)
+ {
+ result = rlcsa.count(rows[i]);
+ occurrences = length(result);
+ #pragma omp critical
+ {
+ if(!isEmpty(result))
+ {
+ found++;
+ if(occurrences <= MAX_OCCURRENCES) { total += occurrences; }
+ else { ignored++; }
+ }
+ total_size += rows[i].length();
+ }
+ if(!isEmpty(result) && occurrences <= MAX_OCCURRENCES)
+ {
+ matches = rlcsa.locate(result);
+ delete[] matches;
+ }
+ }
+ }
+ #else
+ for(i = 0; i < n; i++)
+ {
+ result = rlcsa.count(rows[i]);
+ total_size += rows[i].length();
+ if(!isEmpty(result)) { found++; total += length(result); }
+ matches = rlcsa.locate(result);
+ if(matches) { delete[] matches; }
+ }
+ #endif
+
+ double time = readTimer() - start;
+ std::cout << "Patterns: " << n << " (" << (n / time) << " / sec)" << std::endl;
+ std::cout << "Total size: " << total_size << " bytes (" << (total_size / time) << " / sec)" << std::endl;
+ std::cout << "Matches: " << total << " (" << (total / time) << " / sec)" << std::endl;
+ std::cout << "Found: " << found << std::endl;
+ std::cout << "Time: " << time << " seconds" << std::endl;
+ std::cout << std::endl;
+ std::cout << "Ignored " << ignored << " patterns with more than " << MAX_OCCURRENCES << " occurrences." << std::endl;
+ std::cout << std::endl;
+
+ return 0;
+}
--- /dev/null
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+
+using namespace CSA;
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "LCP sampler" << std::endl;
+ if(argc < 3)
+ {
+ std::cout << "Usage: sample_plcp base_name sample_rate" << std::endl;
+ return 1;
+ }
+
+ std::string base_name = argv[1];
+ std::string lcp_name = base_name + LCP_SAMPLES_EXTENSION;
+ std::cout << "LCP samples: " << lcp_name << std::endl;
+ std::ofstream lcp_file(lcp_name.c_str(), std::ios_base::binary);
+ if(!lcp_file)
+ {
+ std::cerr << "Error creating LCP sample file!" << std::endl;
+ return 2;
+ }
+
+ usint sample_rate = atoi(argv[2]);
+ std::cout << "Sample rate: " << sample_rate << std::endl;
+ std::cout << std::endl;
+ RLCSA* rlcsa = new RLCSA(base_name);
+
+ double start = readTimer();
+ pair_type* sampled_values = 0;
+ usint samples = rlcsa->sampleLCP(sample_rate, sampled_values, true);
+ usint data_size = rlcsa->getSize();
+ delete rlcsa; // Saves memory.
+ LCPSamples lcp(sampled_values, data_size, samples, true, true);
+ lcp.writeTo(lcp_file);
+
+ double time = readTimer() - start;
+ double megabytes = data_size / (double)MEGABYTE;
+ double size = lcp.reportSize() / (double)MEGABYTE;
+ std::cout << megabytes << " megabytes in " << time << " seconds (" << (megabytes / time) << " MB/s)" << std::endl;
+ std::cout << "Sampled LCP size: " << size << " MB" << std::endl;
+ std::cout << std::endl;
+
+ lcp_file.close();
+
+ return 0;
+}
#include <algorithm>
#include <fstream>
#include <iostream>
+#include <vector>
#include "sasamples.h"
#include "misc/utils.h"
+#ifdef MULTITHREAD_SUPPORT
+#ifndef _GLIBCXX_PARALLEL
+#include <mcstl.h>
+#endif
+#endif
+
namespace CSA
{
this->items = indexes->getNumberOfItems();
this->integer_bits = length(this->items - 1);
- this->samples = new FastBitBuffer(sample_file, this->items, this->integer_bits);
+ this->samples = new ReadBuffer(sample_file, this->items, this->integer_bits);
this->buildInverseSamples();
}
-SASamples::SASamples(usint* array, usint data_size, usint sample_rate) :
+SASamples::SASamples(uint* array, usint data_size, usint sample_rate) :
rate(sample_rate),
size(data_size)
{
this->items = (this->size - 1) / this->rate + 1;
this->integer_bits = length(this->items - 1);
- this->samples = new FastBitBuffer(this->items, this->integer_bits);
- DeltaEncoder encoder(SASamples::BLOCK_SIZE);
+ WriteBuffer sample_buffer(this->items, this->integer_bits);
+ DeltaEncoder encoder(INDEX_BLOCK_SIZE);
for(usint i = 0; i < this->size; i++)
{
- if(array[i] % rate == 0)
+ if(array[i] % sample_rate == 0)
{
encoder.setBit(i);
- this->samples->writeItem(array[i] / sample_rate);
+ sample_buffer.writeItem(array[i] / sample_rate);
}
}
+
this->indexes = new DeltaVector(encoder, this->size);
+ this->samples = sample_buffer.getReadBuffer();
+ this->buildInverseSamples();
+}
+
+SASamples::SASamples(uint* inverse, DeltaVector* end_points, usint data_size, usint sample_rate) :
+ rate(sample_rate),
+ items(0)
+{
+ usint sequences = end_points->getNumberOfItems();
+
+ DeltaVector::Iterator iter(*(end_points));
+
+ // Determine the samples, insert them into a vector, and sort them.
+ usint start = 0, end = iter.select(0); // Closed range in padded collection.
+ usint seq_start = 0, seq_end = end; // Closed range in inverse SA.
+ std::vector<pair_type>* vec = new std::vector<pair_type>();
+ for(usint i = 0; i < sequences; i++)
+ {
+ for(usint j = seq_start; j <= seq_end; j += this->rate)
+ {
+ vec->push_back(pair_type(inverse[j] - sequences - 1, this->items));
+ this->items++;
+ }
+ start = nextMultipleOf(this->rate, end);
+ end = iter.selectNext();
+ seq_start = seq_end + 2;
+ seq_end = seq_start + end - start;
+ }
+ #ifdef MULTITHREAD_SUPPORT
+ #ifdef _GLIBCXX_PARALLEL
+ std::sort(vec->begin(), vec->end(), __gnu_parallel::sequential_tag());
+ #else
+ std::sort(vec->begin(), vec->end(), mcstl::sequential_tag());
+ #endif
+ #else
+ std::sort(vec->begin(), vec->end());
+ #endif
+
+ // Compress the samples.
+ this->integer_bits = length(this->items - 1);
+ this->size = end + 1;
+ WriteBuffer sample_buffer(this->items, this->integer_bits);
+ DeltaEncoder encoder(INDEX_BLOCK_SIZE);
+ for(usint i = 0; i < this->items; i++)
+ {
+ pair_type sample = (*vec)[i];
+ encoder.setBit(sample.first);
+ sample_buffer.writeItem(sample.second);
+ }
+ delete vec;
+ this->indexes = new DeltaVector(encoder, this->size);
+ this->samples = sample_buffer.getReadBuffer();
this->buildInverseSamples();
}
-SASamples::SASamples(SASamples& index, SASamples& increment, usint* positions, usint number_of_positions) :
+SASamples::SASamples(SASamples& index, SASamples& increment, usint* positions, usint number_of_positions, usint number_of_sequences) :
rate(index.rate),
size(index.size + increment.size),
items(index.items + increment.items)
{
- this->mergeSamples(index, increment, positions, number_of_positions);
+ this->mergeSamples(index, increment, positions, number_of_positions, number_of_sequences);
index.indexes = 0;
index.samples = 0;
delete this->inverse_samples;
}
-void
-SASamples::writeTo(std::ofstream& sample_file)
-{
- this->indexes->writeTo(sample_file);
- this->samples->writeBuffer(sample_file, false);
-}
-
//--------------------------------------------------------------------------
-usint
-SASamples::inverseSA(usint value)
-{
- if(value >= this->size) { return this->size; }
- usint i = this->inverse_samples->readItem(value / this->rate);
- return this->indexes->select(i);
-}
-
-pair_type
-SASamples::getFirstSampleAfter(usint index)
+void
+SASamples::writeTo(std::ofstream& sample_file) const
{
- if(index >= this->size) { return pair_type(this->size, this->size); }
-
- return this->indexes->valueAfter(index);
+ this->indexes->writeTo(sample_file);
+ this->samples->writeBuffer(sample_file);
}
-//--------------------------------------------------------------------------
-
usint
-SASamples::reportSize()
+SASamples::reportSize() const
{
- usint bytes = sizeof(*this) - sizeof(this->indexes);
+ usint bytes = sizeof(*this);
bytes += this->indexes->reportSize();
bytes += this->samples->reportSize();
bytes += this->inverse_samples->reportSize();
void
SASamples::buildInverseSamples()
{
- this->inverse_samples = new FastBitBuffer(this->items, this->integer_bits);
+ WriteBuffer inverse_buffer(this->items, this->integer_bits);
this->samples->goToItem(0);
for(usint i = 0; i < this->items; i++)
{
- this->inverse_samples->goToItem(this->samples->readItem());
- this->inverse_samples->writeItem(i);
+ inverse_buffer.goToItem(this->samples->readItem());
+ inverse_buffer.writeItem(i);
}
+
+ this->inverse_samples = inverse_buffer.getReadBuffer();
}
//--------------------------------------------------------------------------
-
void
-SASamples::mergeSamples(SASamples& index, SASamples& increment, usint* positions, usint n)
+SASamples::mergeSamples(SASamples& index, SASamples& increment, usint* positions, usint n, usint skip)
{
- DeltaVector* first = index.indexes;
- DeltaVector* second = increment.indexes;
- FastBitBuffer* first_samples = index.samples;
- FastBitBuffer* second_samples = increment.samples;
+ DeltaVector::Iterator first(*(index.indexes));
+ DeltaVector::Iterator second(*(increment.indexes));
+ ReadBuffer* first_samples = index.samples;
+ ReadBuffer* second_samples = increment.samples;
- usint first_bit = first->select(0);
+ usint first_bit = first.select(0);
bool first_finished = false;
- usint second_bit = second->select(0);
+ usint second_bit = second.select(0);
usint sum = index.items;
first_samples->goToItem(0);
second_samples->goToItem(0);
- DeltaEncoder encoder(SASamples::BLOCK_SIZE);
+ DeltaEncoder encoder(INDEX_BLOCK_SIZE);
this->integer_bits = length(this->items - 1);
- this->samples = new FastBitBuffer(this->items, this->integer_bits);
+ WriteBuffer sample_buffer(this->items, this->integer_bits);
for(usint i = 0; i < n; i++, first_bit++)
{
- while(!first_finished && first_bit < positions[i])
+ while(!first_finished && first_bit < positions[i] - skip)
{
encoder.setBit(first_bit);
- this->samples->writeItem(first_samples->readItem());
- if(first->hasNext())
+ sample_buffer.writeItem(first_samples->readItem());
+ if(first.hasNext())
{
- first_bit = first->selectNext() + i;
+ first_bit = first.selectNext() + i;
}
else
{
if(i == second_bit) // positions[i] is one
{
- encoder.setBit(positions[i]);
- this->samples->writeItem(second_samples->readItem() + sum);
- second_bit = second->selectNext();
+ encoder.setBit(positions[i] - skip);
+ sample_buffer.writeItem(second_samples->readItem() + sum);
+ second_bit = second.selectNext();
}
}
while(!first_finished)
{
encoder.setBit(first_bit);
- this->samples->writeItem(first_samples->readItem());
- if(!first->hasNext()) { break; }
- first_bit = first->selectNext() + n;
+ sample_buffer.writeItem(first_samples->readItem());
+ if(!first.hasNext()) { break; }
+ first_bit = first.selectNext() + n;
}
delete index.indexes;
delete increment.inverse_samples;
this->indexes = new DeltaVector(encoder, size);
+ this->samples = sample_buffer.getReadBuffer();
}
class SASamples
{
public:
- const static usint BLOCK_SIZE = 16;
+ const static usint INDEX_BLOCK_SIZE = 16;
SASamples(std::ifstream& sample_file, usint sample_rate);
- SASamples(usint* array, usint data_size, usint sample_rate);
+
+ // These assume < 2 GB data. Use the second one when there are multiple sequences.
+ SASamples(uint* array, usint data_size, usint sample_rate);
+ SASamples(uint* inverse, DeltaVector* end_points, usint data_size, usint sample_rate);
+
~SASamples();
// Destroys contents of index and increment.
// We assume index and increment have same sample rate.
- SASamples(SASamples& index, SASamples& increment, usint* positions, usint number_of_positions);
+ // positions must not containt the positions of end of sequence markers.
+ // number_of_sequences is subtracted from each position before the value is used.
+ SASamples(SASamples& index, SASamples& increment, usint* positions, usint number_of_positions, usint number_of_sequences);
- void writeTo(std::ofstream& sample_file);
+ void writeTo(std::ofstream& sample_file) const;
// Returns i such that SA[i] = value.
// If SA[i] is not sampled, returns the next sampled value. (Don't try!)
// Value is actual 0-based suffix array value.
// Returns size if value is too large.
- usint inverseSA(usint value);
+ inline usint inverseSA(usint value) const
+ {
+ if(value >= this->size) { return this->size; }
+ DeltaVector::Iterator iter(*(this->indexes));
+ return iter.select(this->inverse_samples->readItemConst(value / this->rate));
+ }
// Returns the value of ith sample in suffix array order.
- inline usint getSample(usint i)
+ inline usint getSample(usint i) const
{
- return std::min(this->samples->readItem(i) * this->rate, this->size - 1);
+ return std::min(this->samples->readItemConst(i) * this->rate, this->size - 1);
}
// Returns (ind, sample number) where ind >= index or (size, ???).
- pair_type getFirstSampleAfter(usint index);
+ inline pair_type getFirstSampleAfter(usint index) const
+ {
+ DeltaVector::Iterator iter(*(this->indexes));
+ return iter.valueAfter(index);
+ }
- inline usint getSampleRate() { return this->rate; }
- inline usint getNumberOfSamples() { return this->items; }
+ inline usint getSampleRate() const { return this->rate; }
+ inline usint getNumberOfSamples() const { return this->items; }
- usint reportSize();
+ usint reportSize() const;
private:
usint integer_bits;
DeltaVector* indexes;
- FastBitBuffer* samples;
- FastBitBuffer* inverse_samples;
+ ReadBuffer* samples;
+ ReadBuffer* inverse_samples;
void buildInverseSamples();
// Note: contents of original samples are deleted.
- void mergeSamples(SASamples& index, SASamples& increment, usint* positions, usint n);
+ void mergeSamples(SASamples& index, SASamples& increment, usint* positions, usint n, usint skip);
+
+ // These are not allowed.
+ SASamples();
+ SASamples(const SASamples&);
+ SASamples& operator = (const SASamples&);
};
--- /dev/null
+#include <algorithm>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+
+#include "rlcsa.h"
+#include "misc/utils.h"
+
+
+using namespace CSA;
+
+
+int main(int argc, char** argv)
+{
+ std::cout << "Markov chain text generator" << std::endl;
+ if(argc < 5)
+ {
+ std::cout << "Usage: text_generator index_name initial_context context_length output" << std::endl;
+ return 1;
+ }
+
+ std::cout << "Index name: " << argv[1] << std::endl;
+ std::cout << "Initial context: " << argv[2] << std::endl;
+
+ usint context_length = atoi(argv[3]);
+ std::cout << "Context length: " << context_length << std::endl;
+
+ std::cout << "Output: " << argv[4] << std::endl;
+ std::ofstream output(argv[4], std::ios_base::binary);
+ if(!output)
+ {
+ std::cerr << "Error creating output file!" << std::endl;
+ return 2;
+ }
+ std::cout << std::endl;
+
+ RLCSA rlcsa(argv[1]);
+ srand((long)readTimer());
+ uchar* buffer = new uchar[MEGABYTE];
+ usint len = rlcsa.generateText(argv[2], context_length, buffer, MEGABYTE);
+ output.write((char*)buffer, len);
+
+ delete[] buffer;
+ output.close();
+ return 0;
+}
--- /dev/null
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+#include "../misc/definitions.h"
+#include "../misc/utils.h"
+
+
+using namespace CSA;
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "Converting lines into a patterns file in Pizza & Chili format" << std::endl;
+ if(argc < 3)
+ {
+ std::cout << "Usage: convert_patterns input output" << std::endl;
+ return 1;
+ }
+
+ std::cout << "Input file: " << argv[1] << std::endl;
+ std::ifstream input_file(argv[1], std::ios_base::binary);
+ if(!input_file)
+ {
+ std::cerr << "Error opening input file!" << std::endl;
+ return 2;
+ }
+
+ std::cout << "Output file: " << argv[2] << std::endl;
+ std::ofstream output_file(argv[2], std::ios_base::binary);
+ if(!output_file)
+ {
+ std::cerr << "Error creating output file!" << std::endl;
+ return 3;
+ }
+ std::cout << std::endl;
+
+ std::vector<std::string> rows;
+ readRows(input_file, rows, true);
+
+ usint patterns = rows.size();
+ usint pattern_length = (*(rows.begin())).length();
+ std::cout << "Number of patterns: " << patterns << std::endl;
+ std::cout << "Pattern length: " << pattern_length << std::endl;
+ std::cout << std::endl;
+
+ output_file << "# number=" << patterns << " length=" << pattern_length << " file=foo forbidden=" << std::endl;
+ for(std::vector<std::string>::iterator iter = rows.begin(); iter != rows.end(); iter++)
+ {
+ output_file.write((*iter).c_str(), pattern_length);
+ }
+
+ output_file.close();
+ input_file.close();
+ return 0;
+}
--- /dev/null
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include "../misc/definitions.h"
+
+
+using namespace CSA;
+
+
+std::string gapChars = " =-_";
+
+void
+concatenateStrain(std::string &text, std::ifstream &file, usint strain)
+{
+ file.clear();
+ file.seekg(0, std::ios::beg);
+
+ std::string line;
+ while(!file.eof())
+ {
+ getline(file, line);
+ if(line.empty() || line[0] == '>') { continue; }
+ char t = toupper(line[strain * 2]);
+ if(gapChars.find(t) == std::string::npos) { text.push_back(t); }
+ }
+}
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "Extracting concatenated text from alignment" << std::endl;
+ if(argc < 3)
+ {
+ std::cout << "Usage: extract_text alignment_file output_file" << std::endl;
+ return 1;
+ }
+
+ std::cout << "Alignment file: " << argv[1] << std::endl;
+ std::ifstream alignment_file(argv[1], std::ios_base::binary);
+ if(!alignment_file)
+ {
+ std::cerr << "Error opening alignment file!" << std::endl;
+ return 2;
+ }
+
+ std::cout << "Output file: " << argv[2] << std::endl;
+ std::ofstream output_file(argv[2], std::ios_base::binary);
+ if(!output_file)
+ {
+ std::cerr << "Error creating output file!" << std::endl;
+ return 3;
+ }
+ std::cout << std::endl;
+
+ std::string line;
+ std::getline(alignment_file, line); // Skipping header line.
+ std::getline(alignment_file, line);
+ usint strains = line.size() / 2 - 1;
+ std::cout << "Number of strains: " << strains << std::endl;
+
+ line.clear();
+ concatenateStrain(line, alignment_file, 0);
+ usint n = line.size();
+ std::cout << "Reference sequence length: " << n << std::endl;
+
+ for(usint i = 1; i <= strains; i++)
+ {
+ concatenateStrain(line, alignment_file, i);
+ }
+ usint N = line.size() + 1;
+ std::cout << "Total length: " << N << std::endl;
+
+ output_file.write(line.c_str(), N - 1);
+ char end = '\0';
+ output_file.write(&end, 1);
+
+ output_file.close();
+ alignment_file.close();
+ return 0;
+}
--- /dev/null
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#include "../misc/definitions.h"
+#include "../misc/utils.h"
+
+
+using namespace CSA;
+
+
+int
+main(int argc, char** argv)
+{
+ std::cout << "Text splitter" << std::endl;
+ if(argc < 3)
+ {
+ std::cout << "Usage: split_text file parts" << std::endl;
+ return 1;
+ }
+
+ std::string base_name = argv[1];
+ std::cout << "Input file: " << base_name << std::endl;
+ std::ifstream input_file(base_name.c_str(), std::ios_base::binary);
+ if(!input_file)
+ {
+ std::cerr << "Error opening input file!" << std::endl;
+ return 2;
+ }
+
+ usint parts = atoi(argv[2]);
+ usint size = fileSize(input_file);
+ usint part_size = (size + parts - 1) / parts;
+ std::cout << "Parts: " << parts << std::endl;
+ std::cout << "Part size: " << part_size << " bytes" << std::endl;
+ if(parts == 0 || part_size == 0)
+ {
+ std::cerr << "Invalid number of parts!" << std::endl;
+ return 3;
+ }
+ std::cout << std::endl;
+
+ char* buffer = new char[part_size + 1];
+ usint len = 0;
+ for(usint n = parts; n > 0; n /= 10) { len++; }
+
+ for(usint i = 1; i <= parts; i++)
+ {
+ usint llen = len;
+ for(usint num = i; num > 0; num /= 10) { llen--; }
+
+ std::ostringstream out;
+ out << base_name << "." << std::string(llen, '0') << i;
+ std::ofstream part_file(out.str().c_str(), std::ios_base::binary);
+ input_file.read(buffer, part_size);
+ usint n = input_file.gcount();
+ buffer[n] = 0;
+ part_file.write(buffer, n + 1);
+ part_file.close();
+ }
+
+ input_file.close();
+ return 0;
+}
--- /dev/null
+#! /usr/bin/env python
+# -*- coding: iso-8859-15 -*-
+
+import getopt, sys
+
+
+def main():
+ if len(sys.argv) < 4:
+ return
+
+ sequences = 0
+ in_sequence = False
+ part_size = 1048576 * int(sys.argv[2])
+ current_file = 1
+ print "Part size:", part_size
+ print
+
+ if sys.argv[1] == "-":
+ infile = sys.stdin
+ else:
+ infile = open(sys.argv[1], "r")
+ partname = "part"
+
+ start_tag = "<" + sys.argv[3] + ">"
+ end_tag = "</" + sys.argv[3] + ">"
+
+ output = open(partname + "." + str(current_file), "wb")
+ print "Writing part", output.name, "..."
+
+ for line in infile:
+ if in_sequence:
+ if line.find(end_tag) >= 0:
+ output.write("\0")
+ in_sequence = False
+ else:
+ output.write(line)
+ else:
+ if line.find(start_tag) >= 0:
+ if output.tell() >= part_size:
+ output.close()
+ current_file += 1
+ output = open(partname + "." + str(current_file), "wb")
+ print "Writing part", output.name, "..."
+ in_sequence = True
+ sequences += 1
+
+ infile.close()
+ output.close()
+ print
+ print "Sequences: ", sequences
+
+if __name__ == "__main__":
+ main()
CC = g++
LIBCDSPATH = ../libcds
CPPFLAGS = -Wall -ansi -g -I$(LIBCDSPATH)/includes/ -O3 -DNDEBUG
+#CPPFLAGS = -Wall -ansi -g -I$(LIBCDSPATH)/includes/
LIBCDSA = $(LIBCDSPATH)/lib/libcds.a
LIBRLCSA = incbwt/rlcsa.a
LIBLZTRIE = lzindex/lztrie.a
dcover_obs = dcover/difference_cover.o
-TextCollection_obs = TextCollection.o TextCollectionBuilder.o FMIndexBuilder.o FMIndex.o Tools.o \
+TextCollection_obs = TextCollection.o TextCollectionBuilder.o FMIndexBuilder.o RLCSABuilder.o FMIndex.o Tools.o \
TextStorage.o ${LIBRLCSA} ${LIBCDSA} ${LIBLZTRIE} ${LIBSWCSA}
TCDebug_obs = bittree.o rbtree.o dynFMI.o
projects / SXSI / TextCollection.git / commitdiff