Added SWCSA
[SXSI/TextCollection.git] / FMIndex.cpp
diff --git a/FMIndex.cpp b/FMIndex.cpp
new file mode 100644 (file)
index 0000000..297604b
--- /dev/null
@@ -0,0 +1,1309 @@
+/******************************************************************************
+ *   Copyright (C) 2006-2008 by Veli Mäkinen and Niko Välimäki                *
+ *                                                                            *
+ *   FMIndex 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.            *
+ *****************************************************************************/
+#include "FMIndex.h"
+
+//#define DEBUG_MEMUSAGE
+#ifdef DEBUG_MEMUSAGE
+#include "HeapProfiler.h" // FIXME remove
+#endif
+
+#include <iostream>
+#include <map>
+#include <set>
+#include <vector>
+#include <utility>
+#include <stdexcept>
+#include <cassert>
+#include <cstring> // For strlen()
+using std::vector;
+using std::pair;
+using std::make_pair;
+using std::map;
+using std::string;
+
+namespace SXSI
+{
+
+// Save file version info
+const uchar FMIndex::versionFlag = 9;
+
+/**
+ * Constructor inits an empty dynamic FM-index.
+ * Samplerate defaults to TEXTCOLLECTION_DEFAULT_SAMPLERATE.
+ */
+FMIndex::FMIndex(uchar * bwt, ulong length, unsigned samplerate_, 
+                     unsigned numberOfTexts_, ulong maxTextLength_, ulong numberOfSamples_, 
+                     CSA::DeltaVector & notIndexed, const string & niText, char tsType)
+    : n(length), samplerate(samplerate_), alphabetrank(0), sampled(0), suffixes(0), 
+      suffixDocId(0), numberOfTexts(numberOfTexts_), maxTextLength(maxTextLength_), Doc(0)
+{
+    makewavelet(bwt); // Deletes bwt!
+    bwt = 0;
+    // Make sampling tables
+    maketables(numberOfSamples_, tsType, notIndexed, niText);
+}
+
+bool FMIndex::EmptyText(DocId k) const
+{
+    assert(k < (DocId)numberOfTexts); 
+    return false; // Empty texts are not indexed
+}
+
+uchar * FMIndex::GetText(DocId k) const
+{
+    assert(k < (DocId)numberOfTexts);
+
+    return textStorage->GetText(k);
+/*    TextPosition i = k;
+
+    string result;
+    // Reserve average string length to avoid reallocs
+    result.reserve(n/numberOfTexts);
+
+    uchar c = alphabetrank->access(i);
+    while (c != '\0')
+    {
+        result.push_back(c);
+        i = C[c]+alphabetrank->rank(c,i)-1;
+        
+        c = alphabetrank->access(i); // "next" char.
+    }
+    
+    // Convert to uchar (FIXME return string?)
+    i = result.size();
+    uchar* res = new uchar[i+1]; 
+    res[i] = '\0';   
+    for (ulong j = 0; j < i; ++j)
+        res[i-j-1] = result[j];
+        return res;*/
+}
+
+/*
+ * Substring queries are supported via the pointer returned by TextStorage::GetText
+uchar* FMIndex::GetText(DocId k, TextPosition i, TextPosition j) const
+{
+    assert(k < (DocId)numberOfTexts);
+    assert(j < (*textLength)[k]);
+    assert(i <= j);
+
+    ulong textRank = 0;
+
+    // Start position of k'th text
+    ulong start = (*textStartPos)[k];
+
+    return Substring(i + start, j-i+1);
+    }*/
+
+
+
+/******************************************************************
+ * Existential queries
+ */
+bool FMIndex::IsPrefix(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return true;
+
+    TextPosition sp = 0, ep = 0;
+    Search(pattern, m, &sp, &ep);
+
+    // Check for end-marker(s) in result interval
+    if (CountEndmarkers(sp, ep))
+        return true;
+    return false;
+}
+
+bool FMIndex::IsPrefix(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return true;
+
+    TextPosition sp = 0, ep = 0;
+    Search(pattern, m, &sp, &ep);
+
+    // Check for end-marker(s) in result interval
+    if (CountEndmarkers(sp, ep, begin, end))
+        return true;
+    return false;
+}
+
+
+bool FMIndex::IsSuffix(uchar const *pattern) const
+{
+    // Here counting is as fast as IsSuffix():
+    if (CountSuffix(pattern) > 0)
+        return true;
+    return false;
+}
+
+bool FMIndex::IsSuffix(uchar const *pattern, DocId begin, DocId end) const
+{
+    // Here counting is as fast as IsSuffix():
+    if (CountSuffix(pattern, begin, end) > 0)
+        return true;
+    return false;
+}
+
+bool FMIndex::IsEqual(uchar const *pattern) const
+{
+    TextPosition m = std::strlen((char *)pattern);
+    if (m == 0)
+        return false; // No empty texts exists
+
+    TextPosition sp = 0, ep = 0;
+    // Match including end-marker
+    Search(pattern, m+1, &sp, &ep);
+
+    // Check for end-marker(s) in result interval
+    if (CountEndmarkers(sp, ep))
+        return true;
+    return false;
+}
+
+bool FMIndex::IsEqual(uchar const *pattern, DocId begin, DocId end) const
+{
+    TextPosition m = std::strlen((char *)pattern);
+    if (m == 0)
+        return false; // No empty texts exists
+
+    TextPosition sp = 0, ep = 0;
+    // Match including end-marker
+    Search(pattern, m+1, &sp, &ep, begin, end);
+
+    // Check for end-marker(s) in result interval
+    if (CountEndmarkers(sp, ep))
+        return true;
+    return false;
+}
+
+bool FMIndex::IsContains(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return true;
+
+    TextPosition sp = 0, ep = 0;
+    // Just check if pattern exists somewhere
+    ulong count = Search(pattern, m, &sp, &ep);
+    if (count > 0)
+        return true;
+    return false;
+}
+
+bool FMIndex::IsContains(uchar const * pattern, DocId begin, DocId end) const
+{
+    // Here counting is as fast as existential querying
+    if (CountContains(pattern, begin, end) > 0)
+        return true; // FIXME No need to filter result set
+    return false;
+}
+
+bool FMIndex::IsLessThan(uchar const * pattern) const
+{
+    if (CountLessThan(pattern) > 0)
+        return true;
+    return false;
+}
+
+bool FMIndex::IsLessThan(uchar const * pattern, DocId begin, DocId end) const
+{
+    if (CountLessThan(pattern, begin, end) > 0)
+        return true;
+    return false;
+}
+
+/******************************************************************
+ * Counting queries
+ */
+ulong FMIndex::Count(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return 0;
+
+    TextPosition sp = 0, ep = 0;
+    unsigned count = (unsigned) Search(pattern, m, &sp, &ep);
+    return count;
+}
+
+unsigned FMIndex::CountPrefix(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return numberOfTexts;
+
+    TextPosition sp = 0, ep = 0;
+    Search(pattern, m, &sp, &ep);
+
+    // Count end-markers in result interval
+    return CountEndmarkers(sp, ep);
+}
+
+unsigned FMIndex::CountPrefix(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return numberOfTexts;
+
+    TextPosition sp = 0, ep = 0;
+    Search(pattern, m, &sp, &ep);
+
+    // Count end-markers in result interval
+    return CountEndmarkers(sp, ep, begin, end);
+}
+
+unsigned FMIndex::CountSuffix(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return numberOfTexts;
+
+    TextPosition sp = 0, ep = 0;
+    // Search with end-marker
+    unsigned count = (unsigned) Search(pattern, m+1, &sp, &ep);
+    return count;
+}
+
+unsigned FMIndex::CountSuffix(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return numberOfTexts;
+
+    TextPosition sp = 0, ep = 0;
+    // Search with end-marker
+    unsigned count = (unsigned) Search(pattern, m+1, &sp, &ep, begin, end);
+    return count;
+}
+
+unsigned FMIndex::CountEqual(uchar const *pattern) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return 0; // No empty texts. 
+
+    TextPosition sp = 0, ep = 0;
+    // Match including end-marker
+    Search(pattern, m+1, &sp, &ep);
+
+    // Count end-markers in result interval
+    return CountEndmarkers(sp, ep);
+}
+
+unsigned FMIndex::CountEqual(uchar const *pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return 0; // No empty texts. 
+
+    TextPosition sp = 0, ep = 0;
+    // Match including end-marker
+    Search(pattern, m+1, &sp, &ep, begin, end);
+
+    // Count end-markers in result interval
+    return CountEndmarkers(sp, ep); // Already within [begin, end]
+}
+
+unsigned FMIndex::CountContains(uchar const * pattern) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return numberOfTexts; // Total number of texts. 
+
+    // Here counting is as slow as fetching the result set
+    // because we have to filter out occ's that fall within same document.
+    TextCollection::document_result result = Contains(pattern);
+    return result.size();
+}
+
+unsigned FMIndex::CountContains(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return numberOfTexts; // Total number of texts. 
+
+    // Here counting is as slow as fetching the result set
+    // because we have to filter out occ's that fall within same document.
+    TextCollection::document_result result = Contains(pattern, begin, end);
+    return result.size();
+}
+
+// Less than or equal
+unsigned FMIndex::CountLessThan(uchar const * pattern) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return 0; // No empty texts. 
+
+    TextPosition sp = 0, ep = 0;
+    SearchLessThan(pattern, m, &sp, &ep);
+
+    // Count end-markers in result interval
+    return CountEndmarkers(sp, ep);
+}
+
+unsigned FMIndex::CountLessThan(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return 0; // No empty texts. 
+
+    TextPosition sp = 0, ep = 0;
+    SearchLessThan(pattern, m, &sp, &ep);
+
+    // Count end-markers in result interval
+    return CountEndmarkers(sp, ep, begin, end);
+}
+
+/**
+ * Document reporting queries
+ */
+TextCollection::document_result FMIndex::Prefix(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // FIXME Should return all 1...k
+
+    TextPosition sp = 0, ep = 0;
+    Search(pattern, m, &sp, &ep);
+    
+    // Iterate through end-markers in [sp,ep]:
+    return EnumerateEndmarkers(sp, ep);
+}
+
+TextCollection::document_result FMIndex::Prefix(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // FIXME Should return all 1...k
+
+    TextPosition sp = 0, ep = 0;
+    Search(pattern, m, &sp, &ep);
+
+    // Return end-markers in [sp,ep] and [begin, end]:
+    return EnumerateEndmarkers(sp, ep, begin, end);
+}
+
+TextCollection::document_result FMIndex::Suffix(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // FIXME Should return all 1...k
+
+    TextPosition sp = 0, ep = 0;
+    // Search with end-marker
+    Search(pattern, m+1, &sp, &ep);
+    TextCollection::document_result result;
+    result.reserve(ep-sp+1); // Try to avoid reallocation.
+
+    // Check each occurrence
+    for (; sp <= ep; ++sp)
+    {
+        TextPosition i = sp;
+
+        uchar c = alphabetrank->access(i);
+        while (c != '\0' && !sampled->access(i))
+        {
+            i = C[c]+alphabetrank->rank(c,i)-1;
+            c = alphabetrank->access(i);
+        }
+        // Assert: c == '\0'  OR  sampled->IsBitSet(i)
+
+        if (c == '\0')
+        {
+            // Rank among the end-markers in BWT
+            unsigned endmarkerRank = alphabetrank->rank(0, i) - 1;
+            result.push_back(Doc->access(endmarkerRank));
+        }
+        else // Sampled position
+        {
+            DocId docId = (*suffixDocId)[sampled->rank1(i)-1];
+            result.push_back(docId);
+        }
+    }
+    
+    return result;
+}
+
+TextCollection::document_result FMIndex::Suffix(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // FIXME Should return all 1...k
+
+    TextPosition sp = 0, ep = 0;
+    // Search with end-marker
+    Search(pattern, m+1, &sp, &ep, begin, end);
+    TextCollection::document_result result;
+    result.reserve(ep-sp+1); // Try to avoid reallocation.
+
+    // Check each occurrence, already within [begin, end]
+    for (; sp <= ep; ++sp)
+    {
+        TextPosition i = sp;
+
+        uchar c = alphabetrank->access(i);
+        while (c != '\0' && !sampled->access(i))
+        {
+            i = C[c]+alphabetrank->rank(c,i)-1;
+            c = alphabetrank->access(i);
+        }
+        // Assert: c == '\0'  OR  sampled->IsBitSet(i)
+
+        if (c == '\0')
+        {
+            // Rank among the end-markers in BWT
+            unsigned endmarkerRank = alphabetrank->rank(0, i) - 1;
+            result.push_back(Doc->access(endmarkerRank));
+        }
+        else // Sampled position
+        {
+            DocId docId = (*suffixDocId)[sampled->rank1(i)-1];
+            result.push_back(docId);
+        }
+    }
+    
+    return result;
+}
+
+
+TextCollection::document_result FMIndex::Equal(uchar const *pattern) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // FIXME Should return all empty texts
+
+    TextPosition sp = 0, ep = 0;
+    // Match including end-marker
+    Search(pattern, m+1, &sp, &ep);
+
+    // Report end-markers in result interval
+    return EnumerateEndmarkers(sp, ep);
+}
+
+TextCollection::document_result FMIndex::Equal(uchar const *pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char const *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // FIXME Should return all empty texts
+
+    TextPosition sp = 0, ep = 0;
+    // Match including end-marker
+    Search(pattern, m+1, &sp, &ep, begin, end);
+
+    // Report end-markers in result interval
+    return EnumerateEndmarkers(sp, ep, begin, end);
+}
+
+
+TextCollection::document_result FMIndex::Contains(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result();
+
+    TextPosition sp = 0, ep = 0;
+    // Search all occurrences
+    Search(pattern, m, &sp, &ep);
+
+    // We want unique document indentifiers, using std::set to collect them
+    std::set<DocId> resultSet;
+    EnumerateDocuments(resultSet, sp, ep);
+    
+    // Convert std::set to std::vector
+    TextCollection::document_result result(resultSet.begin(), resultSet.end());
+    return result;
+}
+
+TextCollection::document_result FMIndex::Contains(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result();
+
+    TextPosition sp = 0, ep = 0;
+    // Search all occurrences
+    Search(pattern, m, &sp, &ep);
+
+    // We want unique document indentifiers, using std::set to collect them
+    std::set<DocId> resultSet;
+    EnumerateDocuments(resultSet, sp, ep, begin, end);
+    
+    // Convert std::set to std::vector
+    TextCollection::document_result result(resultSet.begin(), resultSet.end());
+    return result;
+}
+
+
+/**
+ ***
+* * FIXME Lessthan or equal
+ */
+TextCollection::document_result FMIndex::LessThan(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // empty result set
+
+    TextPosition sp = 0, ep = 0;
+    SearchLessThan(pattern, m, &sp, &ep);
+
+    // Report end-markers in result interval
+    return EnumerateEndmarkers(sp, ep);
+}
+
+TextCollection::document_result FMIndex::LessThan(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // empty result set
+
+    TextPosition sp = 0, ep = 0;
+    SearchLessThan(pattern, m, &sp, &ep);
+
+    // Iterate through end-markers in [sp,ep] and [begin, end]:
+    return EnumerateEndmarkers(sp, ep, begin, end);
+}
+
+
+TextCollection::document_result FMIndex::KMismaches(uchar const * pattern, unsigned k) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // empty result set
+
+    suffix_range_vector ranges;
+    kmismatches(ranges, pattern, 0, n-1, m, k);
+    std::set<DocId> resultSet;
+
+    for (suffix_range_vector::iterator it = ranges.begin(); it != ranges.end(); ++it)
+        // Iterate through docs in [sp,ep]:
+        EnumerateDocuments(resultSet, (*it).first, (*it).second);
+
+    // Convert std::set to std::vector
+    TextCollection::document_result result(resultSet.begin(), resultSet.end());
+    return result;
+}
+
+TextCollection::document_result FMIndex::KErrors(uchar const * pattern, unsigned k) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::document_result(); // empty result set
+
+    suffix_range_vector ranges;
+    ulong *dd = new ulong[m+1];
+    for (ulong i=0;i<m+1;i++)
+        dd[i]=i;
+    kerrors(ranges, pattern, 0, n-1, m+k, k, dd, m);
+    delete [] dd;
+    
+    std::set<DocId> resultSet;
+    for (suffix_range_vector::iterator it = ranges.begin(); it != ranges.end(); ++it)
+        // Iterate through docs in [sp,ep]:
+        EnumerateDocuments(resultSet, (*it).first, (*it).second);
+
+    // Convert std::set to std::vector
+    TextCollection::document_result result(resultSet.begin(), resultSet.end());
+    return result;
+}
+
+
+/**
+ * Full result set queries
+ */
+TextCollection::full_result FMIndex::FullContains(uchar const * pattern) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return full_result(); // FIXME Throw exception?
+
+    TextPosition sp = 0, ep = 0;
+    // Search all occurrences
+    Search(pattern, m, &sp, &ep);
+    full_result result;
+    result.reserve(ep-sp+1); // Try to avoid reallocation.
+    EnumeratePositions(result, sp, ep);
+    
+    return result;
+}
+
+TextCollection::full_result FMIndex::FullContains(uchar const * pattern, DocId begin, DocId end) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return full_result(); // FIXME Throw exception?
+
+    TextPosition sp = 0, ep = 0;
+    // Search all occurrences
+    Search(pattern, m, &sp, &ep);
+    full_result result;
+    result.reserve(ep-sp+1); // Try to avoid reallocation.
+    EnumeratePositions(result, sp, ep, begin, end);
+    
+    return result;
+}
+
+TextCollection::full_result FMIndex::FullKMismatches(uchar const * pattern, unsigned k) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::full_result(); // empty result set
+
+    suffix_range_vector ranges;
+    ulong count = kmismatches(ranges, pattern, 0, n-1, m, k);
+
+    TextCollection::full_result result;
+    result.reserve(count); // avoid reallocation.
+    for (suffix_range_vector::iterator it = ranges.begin(); it != ranges.end(); ++it)
+        // Iterate through docs in [sp,ep]:
+        EnumeratePositions(result, (*it).first, (*it).second);
+    return result;
+}
+
+TextCollection::full_result FMIndex::FullKErrors(uchar const * pattern, unsigned k) const
+{
+    TextPosition m = strlen((char *)pattern);
+    if (m == 0)
+        return TextCollection::full_result(); // empty result set
+
+    suffix_range_vector ranges;
+    ulong *dd = new ulong[m+1];
+    for (unsigned i=0;i<m+1;i++)
+        dd[i]=i;
+    ulong count = kerrors(ranges, pattern, 0, n-1, m+k, k, dd, m);
+    delete [] dd;
+
+    TextCollection::full_result result;
+    result.reserve(count); // avoid reallocation.
+    for (suffix_range_vector::iterator it = ranges.begin(); it != ranges.end(); ++it)
+        // Iterate through docs in [sp,ep]:
+        EnumeratePositions(result, (*it).first, (*it).second);
+    return result;
+}
+
+
+/**
+ * Save index to a file handle
+ *
+ * Throws a std::runtime_error exception on i/o error.
+ * First byte that is saved represents the version number of the save file.
+ * In this version files, the data fields are:
+ *  uchar type;
+    uchar versionFlag;
+    TextPosition n;
+    unsigned samplerate;
+    unsigned C[256];
+    TextPosition bwtEndPos;
+    static_sequence * alphabetrank;
+    BSGAP *sampled; 
+    BlockArray *suffixes;
+    BlockArray *suffixDocId;
+    unsigned numberOfTexts;
+    ulong maxTextLength;
+    static_sequence *docId;
+ *
+ * The second parameter (filename's prefix) is ignored.
+ */
+void FMIndex::Save(FILE *file, char const *ignored) const
+{
+    const char type = 'F';
+    // Saving type info:
+    if (std::fwrite(&type, 1, 1, file) != 1)
+        throw std::runtime_error("FMIndex::Save(): file write error (type flag).");
+
+    // Saving version info:
+    if (std::fwrite(&versionFlag, 1, 1, file) != 1)
+        throw std::runtime_error("FMIndex::Save(): file write error (version flag).");
+
+    if (std::fwrite(&(this->n), sizeof(TextPosition), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Save(): file write error (n).");
+    if (std::fwrite(&(this->samplerate), sizeof(unsigned), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Save(): file write error (samplerate).");
+
+    for(ulong i = 0; i < 256; ++i)
+        if (std::fwrite(this->C + i, sizeof(unsigned), 1, file) != 1)
+            throw std::runtime_error("FMIndex::Save(): file write error (C table).");
+
+    if (std::fwrite(&(this->bwtEndPos), sizeof(TextPosition), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Save(): file write error (bwt end position).");
+    
+    alphabetrank->save(file);
+    sampled->save(file);
+    suffixes->Save(file);
+    suffixDocId->Save(file);
+    
+    if (std::fwrite(&(this->numberOfTexts), sizeof(unsigned), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Save(): file write error (numberOfTexts).");
+    if (std::fwrite(&(this->maxTextLength), sizeof(ulong), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Save(): file write error (maxTextLength).");
+
+    Doc->save(file);
+    textStorage->Save(file);
+    fflush(file);
+}
+
+
+/**
+ * Load index from a file handle
+ *
+ * Throws a std::runtime_error exception on i/o error.
+ * For more info, see FMIndex::Save().
+ * 
+ * index_mode_t is defined in TextCollection.h and
+ * defaults to both the index and "naive" text.
+ * 
+ * Note: Samplerate can not be changed during load.
+ */
+FMIndex::FMIndex(FILE *file, index_mode_t im, unsigned samplerate_)
+    : n(0), samplerate(samplerate_), alphabetrank(0), sampled(0), suffixes(0), 
+      suffixDocId(0), numberOfTexts(0), maxTextLength(0), Doc(0)
+{
+    // NB: Type byte has already been read from input
+
+    uchar verFlag = 0;
+    if (std::fread(&verFlag, 1, 1, file) != 1)
+        throw std::runtime_error("FMIndex::Load(): file read error (version flag).");
+    if (verFlag != FMIndex::versionFlag)
+        throw std::runtime_error("FMIndex::Load(): invalid save file version.");
+
+    if (std::fread(&(this->n), sizeof(TextPosition), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Load(): file read error (n).");
+    if (std::fread(&samplerate, sizeof(unsigned), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Load(): file read error (samplerate).");
+// FIXME samplerate can not be changed during load.
+//    if (this->samplerate == 0)
+//        this->samplerate = samplerate;
+
+    for(ulong i = 0; i < 256; ++i)
+        if (std::fread(this->C + i, sizeof(unsigned), 1, file) != 1)
+            throw std::runtime_error("FMIndex::Load(): file read error (C table).");
+
+    if (std::fread(&(this->bwtEndPos), sizeof(TextPosition), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Load(): file read error (bwt end position).");
+
+    alphabetrank = static_sequence::load(file);
+    if (im == index_mode_text_only) { delete alphabetrank; alphabetrank = 0; }
+
+    sampled = static_bitsequence::load(file);
+    if (im == index_mode_text_only) { delete sampled; sampled = 0; }
+    suffixes = new BlockArray(file);
+    if (im == index_mode_text_only) { delete suffixes; suffixes = 0; }
+    suffixDocId = new BlockArray(file);
+    if (im == index_mode_text_only) { delete suffixDocId; suffixDocId = 0; }
+
+    if (std::fread(&(this->numberOfTexts), sizeof(unsigned), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Load(): file read error (numberOfTexts).");
+    if (std::fread(&(this->maxTextLength), sizeof(ulong), 1, file) != 1)
+        throw std::runtime_error("FMIndex::Load(): file read error (maxTextLength).");
+
+    Doc = new ArrayDoc(file); //static_sequence::load(file);
+    if (im == index_mode_text_only) { delete Doc; Doc = 0; }
+
+    textStorage = TextStorage::Load(file);
+
+    // FIXME Construct data structures with new samplerate
+    //maketables(); 
+}
+
+
+
+/**
+ * Rest of the functions follow...
+ */
+ulong FMIndex::searchPrefix(uchar const *pattern, ulong i, ulong *sp, ulong *ep) const
+{
+    int c;
+    while (*sp<=*ep && i>=1) 
+    {
+        c = (int)pattern[--i];
+        *sp = C[c]+alphabetrank->rank(c,*sp-1);
+        *ep = C[c]+alphabetrank->rank(c,*ep)-1;
+    }
+    if (*sp<=*ep)
+        return *ep - *sp + 1;
+    else
+        return 0;
+}
+
+
+ulong FMIndex::kmismatches(suffix_range_vector &result, uchar const *pattern, ulong sp, ulong ep, ulong j, unsigned k) const
+{
+    if (sp>ep) return 0;
+    if (j == 0)
+    {
+        result.push_back(std::make_pair(sp,ep));
+        return ep-sp+1;
+    }
+    int c;
+    ulong spnew;
+    ulong epnew;    
+    int knew;
+    ulong sum=0;
+    if (k==0)
+    {
+        sum = searchPrefix(pattern, j, &sp, &ep);
+        if (sp<=ep)
+            result.push_back(std::make_pair(sp, ep));
+        return sum;
+    }
+    vector<int> chars = alphabetrank->accessAll(sp, ep);
+    for (vector<int>::iterator it = chars.begin(); it != chars.end(); ++it) 
+    {
+        if (*it == 0)
+            continue; // skip '\0'
+        c = *it;
+        spnew = C[c]+alphabetrank->rank(c,sp-1);
+        epnew = C[c]+alphabetrank->rank(c,ep)-1;
+        if (c!=pattern[j-1]) knew = (int)k-1; else knew = k;
+        if (knew>=0) sum += kmismatches(result, pattern, spnew, epnew, j-1, knew);         
+    }
+    return sum;       
+}
+
+//first call kerrors(pattern,1,n,m+k,k,d,m), where d[i]=i
+ulong FMIndex::kerrors(suffix_range_vector &result, uchar const *pattern, ulong sp, ulong ep, ulong j, unsigned k, ulong const *d, ulong m) const
+{
+    ulong sum=0;
+    if (d[m]<=k) // range of suffixes with at most k-errors found
+    {
+        if (sp<=ep)
+            result.push_back(std::make_pair(sp, ep));
+        sum += (sp<=ep)?ep-sp+1:0ul;
+    }
+    if (sp>ep || j==0) 
+        return sum;
+    ulong *dnew = new ulong[m+1];       
+    int c;
+    ulong spnew;
+    ulong p,lowerbound;
+    ulong epnew;    
+    vector<int> chars = alphabetrank->accessAll(sp, ep);
+    for (vector<int>::iterator it = chars.begin(); it != chars.end(); ++it) 
+    { 
+        if (*it == 0)
+            continue; // skip '\0'
+        c = *it;
+        spnew = C[c]+alphabetrank->rank(c,sp-1);
+        epnew = C[c]+alphabetrank->rank(c,ep)-1;
+        if (spnew>epnew) continue;
+        dnew[0]=m+k-j+1;
+        lowerbound=k+1;
+        for (p=1; p<=m; p++) {
+            dnew[p]=myminofthree(d[p]+1,dnew[p-1]+1,(c==pattern[m-p])?d[p-1]:(d[p-1]+1));
+            if (dnew[p]<lowerbound)
+                lowerbound = dnew[p];
+        }
+        if (lowerbound<=k)
+            sum += kerrors(result, pattern, spnew, epnew, j-1, k,dnew,m);         
+    }
+    delete [] dnew;
+    return sum;       
+}
+
+
+ulong FMIndex::Search(uchar const * pattern, TextPosition m, TextPosition *spResult, TextPosition *epResult) const 
+{
+    // use the FM-search replacing function Occ(c,1,i) with alphabetrank->rank(c,i)
+    int c = (int)pattern[m-1]; 
+    TextPosition i=m-1;
+    TextPosition sp = C[c];
+    TextPosition ep = C[c+1]-1;
+    while (sp<=ep && i>=1) 
+    {
+//         printf("i = %lu, c = %c, sp = %lu, ep = %lu\n", i, pattern[i], sp, ep);
+        c = (int)pattern[--i];
+        sp = C[c]+alphabetrank->rank(c,sp-1);
+        ep = C[c]+alphabetrank->rank(c,ep)-1;
+    }
+    *spResult = sp;
+    *epResult = ep;
+    if (sp<=ep)
+        return ep - sp + 1;
+    else
+        return 0;
+}
+
+ulong FMIndex::Search(uchar const * pattern, TextPosition m, TextPosition *spResult, TextPosition *epResult, DocId begin, DocId end) const 
+{
+    // use the FM-search replacing function Occ(c,1,i) with alphabetrank->rank(c,i)
+    int c = (int)pattern[m-1]; 
+    assert(c == 0); // Start from endmarkers
+    TextPosition i=m-1;
+    TextPosition sp = begin;
+    TextPosition ep = end;
+    while (sp<=ep && i>=1) 
+    {
+//         printf("i = %lu, c = %c, sp = %lu, ep = %lu\n", i, pattern[i], sp, ep);
+        c = (int)pattern[--i];
+        sp = C[c]+alphabetrank->rank(c,sp-1);
+        ep = C[c]+alphabetrank->rank(c,ep)-1;
+    }
+    *spResult = sp;
+    *epResult = ep;
+    if (sp<=ep)
+        return ep - sp + 1;
+    else
+        return 0;
+}
+
+
+ulong FMIndex::SearchLessThan(uchar const * pattern, TextPosition m, TextPosition *spResult, TextPosition *epResult) const 
+{
+    // use the FM-search replacing function Occ(c,1,i) with alphabetrank->rank(c,i)
+    uint c = (int)pattern[m-1]; 
+    TextPosition i=m-1;
+    TextPosition sp = 1;
+    TextPosition ep = C[c+1]-1;
+    while (sp<=ep && i>=1) 
+    {
+//         printf("i = %lu, c = %c, sp = %lu, ep = %lu\n", i, pattern[i], sp, ep);
+        c = (int)pattern[--i];
+        uint result = alphabetrank->rankLessThan(c,ep);
+        if (result == ~0u)
+            ep = 0;
+        else
+            ep = C[c]+result-1;
+    }
+    *spResult = sp;
+    *epResult = ep;
+    if (sp<=ep)
+        return ep - sp + 1;
+    else
+        return 0;
+}
+
+
+FMIndex::~FMIndex() {
+    delete alphabetrank;       
+    delete sampled;
+    delete suffixes;
+    delete suffixDocId;
+    delete Doc;
+    delete textStorage;
+}
+
+void FMIndex::makewavelet(uchar *bwt)
+{
+    ulong i, min = 0,
+             max;
+    for (i=0;i<256;i++)
+        C[i]=0;
+    for (i=0;i<n;++i)
+        C[(int)bwt[i]]++;
+    for (i=0;i<256;i++)
+        if (C[i]>0) {min = i; break;}          
+    for (i=255;i>=min;--i)
+        if (C[i]>0) {max = i; break;}
+    
+    ulong prev=C[0], temp;
+    C[0]=0;
+    for (i=1;i<256;i++) {          
+        temp = C[i];
+        C[i]=C[i-1]+prev;
+        prev = temp;
+    }
+//    this->codetable = node::makecodetable(bwt,n);
+//    alphabetrank = new THuffAlphabetRank(bwt,n, this->codetable,0);   
+//    delete [] bwt;
+    //alphabetrank = new RLWaveletTree(bwt, n); // Deletes bwt!
+//  std::cerr << "heap usage: " << HeapProfiler::GetHeapConsumption()/(1024*1024) << " Mbytes" << std::endl;
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "max heap usage before WT: " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes" << std::endl;
+    HeapProfiler::ResetMaxHeapConsumption(); 
+#endif
+
+    alphabet_mapper * am = new alphabet_mapper_none();
+    static_bitsequence_builder * bmb = new static_bitsequence_builder_brw32(8); //rrr02(8); // FIXME samplerate?
+    wt_coder * wtc = new wt_coder_huff(bwt,n,am);//binary(bwt,n,am); // FIXME Huffman shape
+    alphabetrank = new static_sequence_wvtree(bwt,n,wtc,bmb,am);
+    delete bmb;
+    bwt = 0; // already deleted
+   
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "heap usage after WT: " << HeapProfiler::GetHeapConsumption()/(1024*1024) << " Mbytes" << std::endl;
+    std::cerr << "max heap usage after WT: " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes" << std::endl;
+#endif
+}
+
+void FMIndex::maketables(ulong sampleLength, char tsType, CSA::DeltaVector & notIndexed, const string & niText)
+{
+    // Calculate BWT end-marker position (of last inserted text)
+    {
+        ulong i = 0;
+        uint alphabetrank_i_tmp = 0;
+        uchar c  = alphabetrank->access(i, alphabetrank_i_tmp);
+        while (c != '\0')
+        {
+            i = C[c]+alphabetrank_i_tmp-1;
+            c = alphabetrank->access(i, alphabetrank_i_tmp);
+        }
+
+        this->bwtEndPos = i;
+    }
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "heap usage before BWT traverse: " << HeapProfiler::GetHeapConsumption()/(1024*1024) << " / " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes, " <<  HeapProfiler::GetHeapConsumption() << " / " <<  HeapProfiler::GetMaxHeapConsumption() << std::endl;
+    HeapProfiler::ResetMaxHeapConsumption();
+#endif
+
+    // Build up array for text starting positions
+//    BlockArray* textStartPos = new BlockArray(numberOfTexts, Tools::CeilLog2(this->n));
+//    (*textStartPos)[0] = 0; 
+
+    // Mapping from end-markers to doc ID's:
+    unsigned logNumberOfTexts = Tools::CeilLog2(numberOfTexts);
+//    uint *endmarkerDocId = new uint[(numberOfTexts * logNumberOfTexts)/(8*sizeof(uint)) + 1];
+    BlockArray *endmarkerDocId = new BlockArray(numberOfTexts, logNumberOfTexts);
+
+    BlockArray* positions = new BlockArray(sampleLength, Tools::CeilLog2(this->n));
+    uint *sampledpositions = new uint[n/(sizeof(uint)*8)+1];
+    for (ulong i = 0; i < n / (sizeof(uint)*8) + 1; i++)
+        sampledpositions[i] = 0;
+    
+    ulong x,p=bwtEndPos;
+    ulong sampleCount = 0;
+    // Keeping track of text position of prev. end-marker seen
+    ulong posOfSuccEndmarker = n-1;
+    DocId textId = numberOfTexts;
+    ulong ulongmax = 0;
+    ulongmax--;
+    uint alphabetrank_i_tmp =0;
+
+    // Text length = n + number of bytes not indexed.
+    TextStorageBuilder tsbuilder(n + niText.length());
+    ulong tsb_i = n + niText.length(); // Iterator from text length to 0.
+    string::const_reverse_iterator nit_i = niText.rbegin(); // Iterator through non-indexed texts
+
+    for (ulong i=n-1;i<ulongmax;i--) {
+        // i substitutes SA->GetPos(i)
+        x=(i==n-1)?0:i+1;
+
+        uchar c = alphabetrank->access(p, alphabetrank_i_tmp);
+
+        tsbuilder[--tsb_i] = c; // Build TextStorage
+
+        if ((posOfSuccEndmarker - i) % samplerate == 0 && c != '\0')
+        {
+            set_field(sampledpositions,1,p,1);
+            (*positions)[sampleCount] = p;
+            sampleCount ++;
+        }
+
+        if (c == '\0')
+        {
+            unsigned prevTextId = textId; // Cache textId value.
+            --textId; 
+            /**
+             * At first c == '\0' it holds that (prevTextId == numberOfTexts), thus,
+             * we have to search for the first text that is actually *indexed*
+             * to get correct prevTextId.
+             */
+            if (prevTextId == numberOfTexts)
+            {
+                prevTextId = 0;
+                while (notIndexed.isSet(prevTextId))
+                    ++ prevTextId;
+                // Now prevTextId points to the first indexed Doc ID.
+            }
+
+            /**
+             * Insert non-indexed texts
+             */
+            while (notIndexed.isSet(textId))
+            {
+                do {
+                    tsbuilder[tsb_i] = *nit_i;
+                    -- tsb_i;
+                    ++ nit_i;
+                } while (nit_i != niText.rend() && *nit_i != '\0');
+                
+                tsbuilder[tsb_i] = '\0';
+
+                if (textId == 0)
+                    break;
+                --textId;
+            }
+            
+            // Record the order of end-markers in BWT:
+            ulong endmarkerRank = alphabetrank_i_tmp - 1;
+            //set_field(endmarkerDocId, logNumberOfTexts, endmarkerRank, (textId + 1) % numberOfTexts);
+            (*endmarkerDocId)[endmarkerRank] = prevTextId % numberOfTexts;
+
+            // Store text length and text start position:
+            if (textId < (DocId)numberOfTexts - 1)
+            {
+//                (*textStartPos)[textId + 1] = x;  // x-1 is text position of end-marker.
+
+                posOfSuccEndmarker = 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); 
+        }
+        else // Now c != '\0', do LF-mapping:
+            p = C[c]+alphabetrank_i_tmp-1;
+    }
+    while (textId > 0 && notIndexed.isSet(textId-1))
+    {
+        do {
+            -- tsb_i;
+            tsbuilder[tsb_i] = *nit_i;
+            ++ nit_i;
+        } while (nit_i != niText.rend() && *nit_i != '\0');
+        --textId;
+    }
+    assert(textId == 0);
+    assert(tsb_i == 0);
+    assert(nit_i == niText.rend());
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "heap usage before tsbuilder init: " << HeapProfiler::GetHeapConsumption()/(1024*1024) << " / " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes, " <<  HeapProfiler::GetHeapConsumption() << " / " <<  HeapProfiler::GetMaxHeapConsumption() << std::endl;
+    HeapProfiler::ResetMaxHeapConsumption();
+#endif
+
+    textStorage = tsbuilder.InitTextStorage(tsType);
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "heap usage after tsbuilder init: " << HeapProfiler::GetHeapConsumption()/(1024*1024) << " / " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes, " <<  HeapProfiler::GetHeapConsumption() << " / " <<  HeapProfiler::GetMaxHeapConsumption() << std::endl;
+    HeapProfiler::ResetMaxHeapConsumption();
+#endif
+
+    sampled = new static_bitsequence_rrr02(sampledpositions, n, 16);
+    delete [] sampledpositions;
+    assert(sampleCount == sampleLength);
+    assert(sampled->rank1(n-1) == sampleLength);
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "heap usage after sampled bit vector: " << HeapProfiler::GetHeapConsumption()/(1024*1024) << " / " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes, " <<  HeapProfiler::GetHeapConsumption() << " / " <<  HeapProfiler::GetMaxHeapConsumption() << std::endl;
+    HeapProfiler::ResetMaxHeapConsumption();
+#endif
+
+    // Suffixes store an offset from the text start position
+    suffixes = new BlockArray(sampleLength, Tools::CeilLog2(maxTextLength));
+    suffixDocId = new BlockArray(sampleLength, Tools::CeilLog2(numberOfTexts));
+
+    x = n + niText.length() - 2;
+    textId = numberOfTexts - 1;
+    posOfSuccEndmarker = x + 1;
+    for(ulong i = 0; i < sampleLength; i ++) {
+        // Find next sampled text position
+        while ((posOfSuccEndmarker - x) % samplerate != 0 
+               || notIndexed.isSet(textId)) // Loop over non-indexed
+        {
+            --x;
+            assert(x != ~0lu);
+            if (textStorage->IsEndmarker(x))
+            {
+                posOfSuccEndmarker = x--;
+                -- textId;
+            }
+        }
+        assert((*positions)[i] < n);
+        ulong j = sampled->rank1((*positions)[i]);
+
+        assert(j != 0); // if (j==0) j=sampleLength;
+        
+        TextPosition textPos = (x==n-1)?0:x+1;
+        (*suffixDocId)[j-1] = textId; // textStorage->DocIdAtTextPos(textPos);
+        assert(textStorage->DocIdAtTextPos(textPos) == textId);
+
+        assert((*suffixDocId)[j-1] < numberOfTexts);
+        // calculate offset from text start:
+        (*suffixes)[j-1] = textPos - textStorage->TextStartPos((*suffixDocId)[j-1]);
+        --x;
+        if (x != ~0lu && textStorage->IsEndmarker(x))
+        {
+            posOfSuccEndmarker = x--;
+            -- textId;
+        }
+    }
+
+    delete positions;
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "heap usage after sampled arrays: " << HeapProfiler::GetHeapConsumption()/(1024*1024) << " / " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes, " <<  HeapProfiler::GetHeapConsumption() << " / " <<  HeapProfiler::GetMaxHeapConsumption() << std::endl;
+    HeapProfiler::ResetMaxHeapConsumption();
+#endif
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "max heap usage before Doc: " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes" << std::endl;
+    HeapProfiler::ResetMaxHeapConsumption();
+#endif
+
+    /*alphabet_mapper * am = new alphabet_mapper_none();
+    static_bitsequence_builder * bmb = new static_bitsequence_builder_rrr02(32); // FIXME samplerate?
+    Doc = new static_sequence_wvtree_noptrs(endmarkerDocId, numberOfTexts, logNumberOfTexts, bmb, am, true);
+    delete bmb;*/
+    // delete [] endmarkerDocId; // already deleted in static_sequence_wvtree_noptrs!
+
+    Doc = new ArrayDoc(endmarkerDocId);
+
+#ifdef DEBUG_MEMUSAGE
+    std::cerr << "max heap usage after Doc: " << HeapProfiler::GetMaxHeapConsumption()/(1024*1024) << " Mbytes" << std::endl;
+#endif
+}
+
+
+/**
+ * Finds document identifier for given text position
+ *
+ * Starting text position of the document is stored into second parameter.
+ * Binary searching on text starting positions. 
+ */
+TextCollection::DocId FMIndex::DocIdAtTextPos(BlockArray* textStartPos, TextPosition i) const
+{
+    assert(i < n);
+
+    DocId a = 0;
+    DocId b = numberOfTexts - 1;
+    while (a < b)
+    {
+        DocId c = a + (b - a)/2;
+        if ((*textStartPos)[c] > i)
+            b = c - 1;
+        else if ((*textStartPos)[c+1] > i)
+            return c;
+        else
+            a = c + 1;
+    }
+
+    assert(a < (DocId)numberOfTexts);
+    assert(i >= (*textStartPos)[a]);
+    assert(i < (a == (DocId)numberOfTexts - 1 ? n : (*textStartPos)[a+1]));
+    return a;
+}
+
+
+} // namespace SXSI
+