CSA.cpp

   1 /******************************************************************************
   2  *   Copyright (C) 2006-2008 by Veli Mäkinen and Niko Välimäki                *
   3  *                                                                            *
   4  *                                                                            *
   5  *   This program is free software; you can redistribute it and/or modify     *
   6  *   it under the terms of the GNU Lesser General Public License as published *
   7  *   by the Free Software Foundation; either version 2 of the License, or     *
   8  *   (at your option) any later version.                                      *
   9  *                                                                            *
  10  *   This program is distributed in the hope that it will be useful,          *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of           *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            *
  13  *   GNU Lesser General Public License for more details.                      *
  14  *                                                                            *
  15  *   You should have received a copy of the GNU Lesser General Public License *
  16  *   along with this program; if not, write to the                            *
  17  *   Free Software Foundation, Inc.,                                          *
  18  *   51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.            *
  19  *****************************************************************************/
  20 #include "CSA.h"
  21 #include "TextCollection.h"
  22 #include <iostream>
  23 #include <queue>
  24 #include <iomanip>
  25 #include <set>
  26 #include <vector>
  27 #include <utility>
  28 #include <stdexcept>
  29 #include <cassert>
  30 #include <cstring> // For strlen()
  31 using std::vector;
  32 using std::pair;
  33 using std::make_pair;
  34 using std::map;
  35
  36 using SXSI::TextCollection;
  37
  38 // Save file version info
  39 const uchar CSA::versionFlag = 2;
  40
  41 ////////////////////////////////////////////////////////////////////////////
  42 // Class CSA::THuffAlphabetRank
  43
  44 CSA::THuffAlphabetRank::THuffAlphabetRank(uchar *s, TextPosition n, TCodeEntry *codetable, unsigned level) {
  45     left = NULL;
  46     right = NULL;
  47     bitrank = NULL;
  48     ch = s[0];
  49     leaf = false;
  50     this->codetable = codetable;
  51
  52     bool *B = new bool[n];
  53     TextPosition sum=0,i;
  54     /*
  55     for (i=0; i< n; i++) {
  56        printf("%c:", (char)((int)s[i]-128));
  57        for (r=0;r<codetable[(int)s[i]].bits;r++)
  58           if (codetable[(int)s[i]].code & (1u <<r))
  59              printf("1");
  60           else printf("0");
  61        printf("\n");
  62     }
  63     printf("\n");
  64     if (level > 100) return;*/
  65     for (i=0;i<n;i++)
  66        if (codetable[(int)s[i]].code & (1u << level)) {
  67           B[i] = true;
  68           sum++;
  69        }
  70        else B[i] = false;
  71     if (sum==0 || sum==n) {
  72         delete [] B;
  73         leaf = true;
  74         return;
  75     }
  76     uchar *sfirst, *ssecond;
  77     //if (n-sum > 0)
  78         sfirst = new uchar[n-sum];
  79     //if (sum > 0)
  80         ssecond = new uchar[sum];
  81     unsigned j=0,k=0;
  82    for (i=0;i<n;i++)
  83         if (B[i]) ssecond[k++] = s[i];
  84         else sfirst[j++] = s[i];
  85     ulong *Binbits = new ulong[n/W+1];
  86     for (i=0;i<n;i++)
  87         Tools::SetField(Binbits,1,i,B[i]);
  88     delete [] B;
  89     bitrank = new BitRank(Binbits,n,true);
  90     //if (j > 0) {
  91         left = new THuffAlphabetRank(sfirst,j,codetable,level+1);
  92         delete [] sfirst;
  93     //}
  94     //if (k>0) {
  95         right = new THuffAlphabetRank(ssecond,k,codetable,level+1);
  96         delete [] ssecond;
  97     //}
  98 }
  99
 100
 101 bool CSA::THuffAlphabetRank::Test(uchar *s, ulong n) {
 102     // testing that the code works correctly
 103     int C[256];
 104     unsigned i,j;
 105     bool correct=true;
 106     for (j=0;j<256;j++)
 107         C[j] = 0;
 108     for (i=0;i<n;i++) {
 109         C[(int)s[i]]++;
 110         if (C[(int)s[i]] != (int)rank((int)s[i],i)) {
 111         correct = false;
 112         printf("%d (%c): %d<>%d\n",i,(int)s[i]-128,C[(int)s[i]],(int)rank((int)s[i],i));
 113         }
 114     }
 115     return correct;
 116 }
 117
 118 CSA::THuffAlphabetRank::~THuffAlphabetRank() {
 119     if (left!=NULL) delete left;
 120     if (right!=NULL) delete right;
 121     if (bitrank!=NULL)
 122         delete bitrank;
 123 }
 124
 125 /**
 126  * Saving data fields:
 127         BitRank *bitrank;
 128         bool leaf;
 129         uchar ch;
 130         left child;
 131         right child;
 132 */
 133 void CSA::THuffAlphabetRank::Save(FILE *file)
 134 {
 135
 136 }
 137
 138 CSA::THuffAlphabetRank::THuffAlphabetRank(FILE *file)
 139 {
 140
 141
 142 }
 143
 144 /////////////////////////////////////////////f///////////////////////////////
 145 // Class CSA
 146
 147 /**
 148  * Constructor inits an empty dynamic FM-index.
 149  * Samplerate defaults to TEXTCOLLECTION_DEFAULT_SAMPLERATE.
 150  */
 151 CSA::CSA(unsigned samplerate)
 152     : n(0), samplerate(0), alphabetrank(0), sampled(0), suffixes(0), suffixDocId(0), positions(0),
 153       codetable(0), dynFMI(0), numberOfTexts(0), numberOfAllTexts(0), maxTextLength(0), endmarkerDocId(0),
 154       textLength(0), textStartPos(0), emptyTextRank(0)
 155 {
 156     this->samplerate = samplerate;
 157
 158     // FIXME TODO : DynFMI needs distribution of characters before hand
 159     // This will create fully balanced wavelet tree for all chars [0, 255].
 160     uchar temp[256];
 161     for (unsigned i = 0; i < 255; ++i)
 162         temp[i] = i+1;
 163     temp[255] = 0;
 164     dynFMI = new DynFMI(temp, 1, 255, false);
 165
 166     /* Debug code: take char distribution from data.txt.
 167     uchar *temp = Tools::GetFileContents("data.txt", 0);
 168     dynFMI = new DynFMI(temp,1790000, strlen((char *)temp), false);
 169     delete [] temp;*/
 170 }
 171
 172 /**
 173  * Insert new text
 174  *
 175  * Given text must include end-marker.
 176  * Text identifiers are numbered starting from 0.
 177  */
 178 void CSA::InsertText(uchar const * text)
 179 {
 180     // Sanity check:
 181     assert(dynFMI != 0);
 182
 183     TextPosition m = std::strlen((char *)text) + 1;
 184     if (m > maxTextLength)
 185         maxTextLength = m; // Store length of the longest text seen so far.
 186
 187     if (m > 1)
 188     {
 189         this->n += m;
 190         this->numberOfTexts ++;
 191         this->numberOfAllTexts ++;
 192         dynFMI->addText(text, m);
 193     }
 194     else
 195     {
 196         emptyTextId.insert(numberOfAllTexts); // FIXME Using too much space here
 197         this->numberOfAllTexts ++;            // Replace with dynamic bitvector
 198     }
 199 }
 200
 201 void CSA::MakeStatic()
 202 {
 203     // Sanity check:
 204     if (dynFMI == 0)
 205         throw std::runtime_error("CSA::MakeStatic(): Data structure is already static (dynFMI == 0).");
 206
 207     // Bitvector of empty texts
 208     {
 209         //std::cout << std::endl << "texts: " << numberOfTexts << ", all texts " << numberOfAllTexts << ", size : " << emptyTextId.size() <<std::endl;
 210         ulong *tempB = new ulong[numberOfAllTexts/W + 1];
 211         for (ulong i = 0; i < numberOfAllTexts/W + 1; ++i)
 212             tempB[i] = 0;
 213         for (std::set<unsigned>::const_iterator it = emptyTextId.begin(); it != emptyTextId.end(); ++it)
 214             Tools::SetField(tempB, 1, (*it), 1);
 215         emptyTextId.clear();
 216         emptyTextRank = new BSGAP(tempB, numberOfAllTexts, true);
 217     }
 218
 219     uchar *bwt = dynFMI->getBWT();
 220     /*printf("123456789012345678901234567890123456789\n");
 221     for (TextPosition i = 0; i < n; i ++)
 222         if (bwt[i] < 50)
 223             printf("%d", (int)bwt[i]);
 224         else
 225             printf("%c", bwt[i]);
 226             printf("\n");*/
 227
 228
 229 /*    for (TextPosition i = 1; i <= n; i ++)
 230     {
 231         printf("LF[%lu, %c] = %lu\n", i, (*dynFMI)[i], dynFMI->LFmapping(i));
 232     }*/
 233
 234     // Sanity check
 235     assert(numberOfTexts == dynFMI->getCollectionSize());
 236
 237     delete dynFMI;
 238     dynFMI = 0;
 239
 240     makewavelet(bwt); // Deletes bwt!
 241     bwt = 0;
 242
 243     // Calculate BWT end-marker position (of last inserted text)
 244     // and the length of the first text (counter l):
 245     ulong i = 0;
 246     ulong l = 1;
 247
 248     uchar c  = alphabetrank->access(i);
 249     while (c != '\0')
 250     {
 251         i = C[c]+alphabetrank->rank(c, i)-1;
 252         ++l;
 253         c = alphabetrank->access(i);
 254     }
 255     bwtEndPos = i;
 256     assert(bwtEndPos < n);
 257     //printf("bwtEndPos = %lu\n", bwtEndPos);
 258
 259     // Build up arrays for text length and starting positions
 260     // FIXME Remove, temp data
 261     textLength = new BlockArray(numberOfTexts, Tools::CeilLog2(maxTextLength));
 262     textStartPos = new BlockArray(numberOfTexts, Tools::CeilLog2(this->n));
 263     (*textLength)[0] = l;
 264     (*textStartPos)[0] = 0; // Rest of the values are updated in CSA::maketables().
 265
 266     // Make sampling tables
 267     maketables();
 268 }
 269
 270 bool CSA::EmptyText(DocId k) const
 271 {
 272     assert(k < (DocId)numberOfTexts);
 273     if (emptyTextRank->IsBitSet(k))
 274         return true;
 275     return false;
 276 }
 277
 278 uchar* CSA::GetText(DocId k) const
 279 {
 280     assert(k < (DocId)numberOfTexts);
 281
 282     ulong textRank = 0;
 283     if (emptyTextRank->IsBitSet(k, &textRank))
 284     {
 285         uchar* result = new uchar[1];
 286         result[0] = 0;
 287         return result;
 288     }
 289     // Map to non-empty text
 290     k -= textRank; //emptyTextRank->rank(k);
 291
 292     TextPosition i = k;
 293
 294     string result;
 295     // Reserve average string length to avoid reallocs
 296     result.reserve(n/numberOfTexts);
 297
 298     uchar c = alphabetrank->access(i);
 299     while (c != '\0')
 300     {
 301         result.push_back(c);
 302         i = C[c]+alphabetrank->rank(c,i)-1;
 303
 304         c = alphabetrank->access(i); // "next" char.
 305     }
 306
 307     // Convert to uchar (FIXME return string?)
 308     i = result.size();
 309     uchar* res = new uchar[i+1];
 310     res[i] = '\0';
 311     for (ulong j = 0; j < i; ++j)
 312         res[i-j-1] = result[j];
 313     return res;
 314 }
 315 /*
 316  * Not supported
 317 uchar* CSA::GetText(DocId k, TextPosition i, TextPosition j) const
 318 {
 319     assert(k < (DocId)numberOfTexts);
 320     assert(j < (*textLength)[k]);
 321     assert(i <= j);
 322
 323     ulong textRank = 0;
 324     if (emptyTextRank->IsBitSet(k, &textRank))
 325     {
 326         uchar* result = new uchar[1];
 327         result[0] = 0;
 328         return result;
 329     }
 330
 331     // Map to non-empty text
 332     k -= textRank; // emptyTextRank->rank(k);
 333
 334     // Start position of k'th text
 335     ulong start = (*textStartPos)[k];
 336
 337     return Substring(i + start, j-i+1);
 338     }*/
 339
 340 /******************************************************************
 341  * Existential queries
 342  */
 343 bool CSA::IsPrefix(uchar const * pattern) const
 344 {
 345     TextPosition m = strlen((char *)pattern);
 346     if (m == 0)
 347         return true;
 348
 349     TextPosition sp = 0, ep = 0;
 350     Search(pattern, m, &sp, &ep);
 351
 352     // Check for end-marker(s) in result interval
 353     if (CountEndmarkers(sp, ep))
 354         return true;
 355     return false;
 356 }
 357
 358 bool CSA::IsSuffix(uchar const *pattern) const
 359 {
 360     // Here counting is as fast as IsSuffix():
 361     if (CountSuffix(pattern) > 0)
 362         return true;
 363     return false;
 364 }
 365
 366 bool CSA::IsEqual(uchar const *pattern) const
 367 {
 368     TextPosition m = std::strlen((char *)pattern);
 369     if (m == 0)
 370     {
 371         if (numberOfAllTexts - numberOfTexts > 0u)
 372             return true; // Empty texts exists
 373         return false; // No empty texts exists
 374     }
 375
 376     TextPosition sp = 0, ep = 0;
 377     // Match including end-marker
 378     Search(pattern, m+1, &sp, &ep);
 379
 380     // Check for end-marker(s) in result interval
 381     if (CountEndmarkers(sp, ep))
 382         return true;
 383     return false;
 384 }
 385
 386 bool CSA::IsContains(uchar const * pattern) const
 387 {
 388     TextPosition m = strlen((char *)pattern);
 389     if (m == 0)
 390         return true;
 391
 392     TextPosition sp = 0, ep = 0;
 393     // Just check if pattern exists somewhere
 394     ulong count = Search(pattern, m, &sp, &ep);
 395
 396     if (count > 0)
 397         return true;
 398     return false;
 399 }
 400
 401 bool CSA::IsLessThan(uchar const*) const
 402 {
 403     // TODO
 404     assert(0);
 405     return false;
 406 }
 407
 408 /******************************************************************
 409  * Counting queries
 410  */
 411 unsigned CSA::CountPrefix(uchar const * pattern) const
 412 {
 413     TextPosition m = strlen((char *)pattern);
 414     if (m == 0)
 415         return numberOfAllTexts;
 416
 417     TextPosition sp = 0, ep = 0;
 418     Search(pattern, m, &sp, &ep);
 419
 420     // Count end-markers in result interval
 421     return CountEndmarkers(sp, ep);
 422 }
 423
 424 unsigned CSA::CountSuffix(uchar const * pattern) const
 425 {
 426     TextPosition m = strlen((char *)pattern);
 427     if (m == 0)
 428         return numberOfAllTexts;
 429
 430     TextPosition sp = 0, ep = 0;
 431     // Search with end-marker
 432     unsigned count = (unsigned) Search(pattern, m+1, &sp, &ep);
 433
 434     return count;
 435 }
 436
 437 unsigned CSA::CountEqual(uchar const *pattern) const
 438 {
 439     TextPosition m = strlen((char const *)pattern);
 440     if (m == 0)
 441         return numberOfAllTexts - numberOfTexts; // Empty texts.
 442
 443     TextPosition sp = 0, ep = 0;
 444     // Match including end-marker
 445     Search(pattern, m+1, &sp, &ep);
 446
 447     // Count end-markers in result interval
 448     return CountEndmarkers(sp, ep);
 449 }
 450
 451 unsigned CSA::CountContains(uchar const * pattern) const
 452 {
 453     TextPosition m = strlen((char const *)pattern);
 454     if (m == 0)
 455         return numberOfAllTexts; // Total number of texts.
 456
 457     // Here counting is as slow as fetching the result set
 458     // because we would have to filter out occ's that fall within same document.
 459     TextCollection::document_result result = Contains(pattern);
 460     return result.size();
 461 }
 462
 463 unsigned CSA::CountLessThan(uchar const * pattern) const
 464 {
 465     // TODO
 466     assert(0);
 467     return 0;
 468 }
 469
 470 /**
 471  * Document reporting queries
 472  */
 473 TextCollection::document_result CSA::Prefix(uchar const * pattern) const
 474 {
 475     TextPosition m = strlen((char *)pattern);
 476     if (m == 0)
 477         return TextCollection::document_result(); // FIXME Should return all 1...k
 478
 479     TextPosition sp = 0, ep = 0;
 480     Search(pattern, m, &sp, &ep);
 481
 482     TextCollection::document_result result;
 483
 484     // Report end-markers in result interval
 485     unsigned resultSize = CountEndmarkers(sp, ep);
 486     if (resultSize == 0)
 487         return result;
 488
 489     result.reserve(resultSize); // Try to avoid reallocation.
 490
 491     // Iterate through end-markers in [sp,ep]:
 492     unsigned i = 0;
 493     if (sp > 0)
 494         i = alphabetrank->rank(0, sp - 1);
 495     while (resultSize)
 496     {
 497         // End-marker we found belongs to the "preceeding" doc in the collection
 498         DocId docId = ((*endmarkerDocId)[i] + 1) % numberOfTexts;
 499         // Map to doc ID:
 500         docId = emptyTextRank->select0(docId+1);
 501         result.push_back(docId);
 502
 503         -- resultSize;
 504         ++ i;
 505     }
 506
 507     return result;
 508 }
 509
 510 TextCollection::document_result CSA::Suffix(uchar const * pattern) const
 511 {
 512     TextPosition m = strlen((char *)pattern);
 513     if (m == 0)
 514         return TextCollection::document_result(); // FIXME Should return all 1...k
 515
 516     TextPosition sp = 0, ep = 0;
 517     // Search with end-marker
 518     Search(pattern, m+1, &sp, &ep);
 519
 520     TextCollection::document_result result;
 521     result.reserve(ep-sp+1); // Try to avoid reallocation.
 522
 523     ulong sampled_rank_i = 0;
 524     // Check each occurrence
 525     for (; sp <= ep; ++sp)
 526     {
 527         TextPosition i = sp;
 528
 529         uchar c = alphabetrank->access(i);
 530         while (c != '\0' && !sampled->IsBitSet(i, &sampled_rank_i))
 531         {
 532             i = C[c]+alphabetrank->rank(c,i)-1;
 533             c = alphabetrank->access(i);
 534         }
 535         // Assert: c == '\0'  OR  sampled->IsBitSet(i)
 536
 537         if (c == '\0')
 538         {
 539             // Rank among the end-markers in BWT
 540             unsigned endmarkerRank = alphabetrank->rank(0, i) - 1;
 541
 542             // End-marker that we found belongs to the "preceeding" doc in collection:
 543             DocId docId = ((*endmarkerDocId)[endmarkerRank] + 1) % numberOfTexts;
 544             // Map to doc ID:
 545             docId = emptyTextRank->select0(docId+1);
 546             result.push_back(docId);
 547         }
 548         else // Sampled position
 549         {
 550             DocId docId = (*suffixDocId)[sampled_rank_i-1]; //sampled->rank(i)-1];
 551             // Map to doc ID:
 552             docId = emptyTextRank->select0(docId+1);
 553             result.push_back(docId);
 554         }
 555     }
 556
 557     return result;
 558 }
 559
 560 TextCollection::document_result CSA::Equal(uchar const *pattern) const
 561 {
 562     TextPosition m = strlen((char const *)pattern);
 563     if (m == 0)
 564         return TextCollection::document_result(); // FIXME Should return all empty texts
 565
 566     TextPosition sp = 0, ep = 0;
 567     // Match including end-marker
 568     Search(pattern, m+1, &sp, &ep);
 569
 570     TextCollection::document_result result;
 571
 572     // Report end-markers in result interval
 573     unsigned resultSize = CountEndmarkers(sp, ep);
 574     if (resultSize == 0)
 575         return result;
 576
 577     result.reserve(resultSize); // Try to avoid reallocation.
 578
 579     unsigned i = 0;
 580     if (sp > 0)
 581         i = alphabetrank->rank(0, sp - 1);
 582     while (resultSize)
 583     {
 584         // End-marker we found belongs to the "previous" doc in the collection
 585         DocId docId = ((*endmarkerDocId)[i] + 1) % numberOfTexts;
 586         // Map to doc ID:
 587         docId = emptyTextRank->select0(docId+1);
 588         result.push_back(docId);
 589
 590         -- resultSize;
 591         ++ i;
 592     }
 593
 594     return result;
 595 }
 596
 597
 598 TextCollection::document_result CSA::Contains(uchar const * pattern) const
 599 {
 600     TextPosition m = strlen((char *)pattern);
 601     if (m == 0)
 602         return TextCollection::document_result();
 603
 604     TextPosition sp = 0, ep = 0;
 605     // Search all occurrences
 606     Search(pattern, m, &sp, &ep);
 607
 608     // We want unique document indentifiers, using std::set to collect them
 609     std::set<DocId> resultSet;
 610
 611     ulong sampled_rank_i = 0;
 612     // Check each occurrence
 613     for (; sp <= ep; ++sp)
 614     {
 615         TextPosition i = sp;
 616         uchar c = alphabetrank->access(i);
 617         while (c != '\0' && !sampled->IsBitSet(i, &sampled_rank_i))
 618         {
 619             i = C[c]+alphabetrank->rank(c,i)-1; // LF-mapping
 620             c = alphabetrank->access(i);
 621         }
 622         if (c == '\0')
 623         {
 624             // Rank among the end-markers in BWT
 625             unsigned endmarkerRank = alphabetrank->rank(0, i) - 1;
 626
 627             // End-marker that we found belongs to the "preceeding" doc in collection:
 628             DocId docId = ((*endmarkerDocId)[endmarkerRank] + 1) % numberOfTexts;
 629             resultSet.insert(docId);
 630         }
 631         else
 632         {
 633             DocId di = (*suffixDocId)[sampled_rank_i-1]; //sampled->rank(i)-1];
 634             assert((unsigned)di < numberOfTexts);
 635             resultSet.insert(di);
 636         }
 637     }
 638
 639     // Convert std::set to std::vector
 640     TextCollection::document_result result(resultSet.begin(), resultSet.end());
 641     // Map to doc ID's
 642     for (document_result::iterator it = result.begin(); it != result.end(); ++it)
 643         *it = emptyTextRank->select0(*it+1);
 644     return result;
 645 }
 646
 647 TextCollection::document_result CSA::LessThan(uchar const * pattern) const
 648 {
 649     // TODO
 650     assert(0);
 651     return document_result();
 652 }
 653
 654 /**
 655  * Full result set queries
 656  */
 657 TextCollection::full_result CSA::FullContains(uchar const * pattern) const
 658 {
 659     TextPosition m = strlen((char *)pattern);
 660     if (m == 0)
 661         return full_result(); // FIXME Throw exception?
 662
 663     TextPosition sp = 0, ep = 0;
 664     // Search all occurrences
 665     Search(pattern, m, &sp, &ep);
 666
 667     full_result result;
 668     result.reserve(ep-sp+1); // Try to avoid reallocation.
 669
 670     ulong sampled_rank_i = 0;
 671     // Report each occurrence
 672     for (; sp <= ep; ++sp)
 673     {
 674         TextPosition i = sp;
 675         TextPosition dist = 0;
 676         uchar c = alphabetrank->access(i);
 677         while (c != '\0' && !sampled->IsBitSet(i, &sampled_rank_i))
 678         {
 679             i = C[c]+alphabetrank->rank(c,i)-1;
 680             c = alphabetrank->access(i);
 681             ++ dist;
 682         }
 683         if (c == '\0')
 684         {
 685             // Rank among the end-markers in BWT
 686             unsigned endmarkerRank = alphabetrank->rank(0, i) - 1;
 687
 688             // End-marker that we found belongs to the "preceeding" doc in collection:
 689             DocId docId = ((*endmarkerDocId)[endmarkerRank] + 1) % numberOfTexts;
 690             // Map to doc ID:
 691             docId = emptyTextRank->select0(docId+1);
 692             result.push_back(make_pair(docId, dist));
 693         }
 694         else
 695         {
 696             TextPosition textPos = (*suffixes)[sampled_rank_i-1]+dist; //sampled->rank(i)-1] + dist;
 697             DocId docId = (*suffixDocId)[sampled_rank_i-1]; //sampled->rank(i)-1];
 698 //            textPos = textPos - (*textStartPos)[docId]; // Offset inside the text
 699
 700             // Map to doc ID:
 701             docId = emptyTextRank->select0(docId+1);
 702             result.push_back(make_pair(docId, textPos));
 703         }
 704     }
 705
 706     return result;
 707 }
 708
 709
 710 /**
 711  * Save index to a file handle
 712  *
 713  * Throws a std::runtime_error exception on i/o error.
 714  * First byte that is saved represents the version number of the save file.
 715  * In version 2 files, the data fields are:
 716  *  uchar versionFlag;
 717     TextPosition n;
 718     unsigned samplerate;
 719     unsigned C[256];
 720     TextPosition bwtEndPos;
 721     static_sequence * alphabetrank;
 722     BSGAP *sampled;
 723     BlockArray *suffixes;
 724     BlockArray *suffixDocId;
 725     unsigned numberOfTexts;
 726     unsigned numberOfAllTexts;
 727     ulong maxTextLength;
 728     BlockArray *endmarkerDocId;
 729     BSGAP *emptyTextRank;
 730  */
 731 void CSA::Save(FILE *file) const
 732 {
 733     // Saving version info:
 734     if (std::fwrite(&versionFlag, 1, 1, file) != 1)
 735         throw std::runtime_error("CSA::Save(): file write error (version flag).");
 736
 737     if (std::fwrite(&(this->n), sizeof(TextPosition), 1, file) != 1)
 738         throw std::runtime_error("CSA::Save(): file write error (n).");
 739     if (std::fwrite(&(this->samplerate), sizeof(unsigned), 1, file) != 1)
 740         throw std::runtime_error("CSA::Save(): file write error (samplerate).");
 741
 742     for(ulong i = 0; i < 256; ++i)
 743         if (std::fwrite(this->C + i, sizeof(unsigned), 1, file) != 1)
 744             throw std::runtime_error("CSA::Save(): file write error (C table).");
 745
 746     if (std::fwrite(&(this->bwtEndPos), sizeof(TextPosition), 1, file) != 1)
 747         throw std::runtime_error("CSA::Save(): file write error (bwt end position).");
 748
 749     alphabetrank->save(file);
 750     sampled->Save(file);
 751     suffixes->Save(file);
 752     suffixDocId->Save(file);
 753
 754     if (std::fwrite(&(this->numberOfTexts), sizeof(unsigned), 1, file) != 1)
 755         throw std::runtime_error("CSA::Save(): file write error (numberOfTexts).");
 756     if (std::fwrite(&(this->numberOfAllTexts), sizeof(unsigned), 1, file) != 1)
 757         throw std::runtime_error("CSA::Save(): file write error (numberOfAllTexts).");
 758     if (std::fwrite(&(this->maxTextLength), sizeof(ulong), 1, file) != 1)
 759         throw std::runtime_error("CSA::Save(): file write error (maxTextLength).");
 760
 761     endmarkerDocId->Save(file);
 762     emptyTextRank->Save(file);
 763 }
 764
 765
 766 /**
 767  * Load index from a file handle
 768  *
 769  * Throws a std::runtime_error exception on i/o error.
 770  * For more info, see CSA::Save().
 771  */
 772 void CSA::Load(FILE *file, unsigned samplerate)
 773 {
 774     // Delete everything
 775     delete dynFMI;       dynFMI = 0;
 776     delete alphabetrank; alphabetrank = 0;
 777     delete sampled;      sampled = 0;
 778     delete suffixes;     suffixes = 0;
 779     delete suffixDocId;  suffixDocId = 0;
 780     delete positions; positions = 0;
 781     delete [] codetable; codetable = 0;
 782
 783     delete endmarkerDocId; endmarkerDocId = 0;
 784     delete emptyTextRank;  emptyTextRank = 0;
 785     // FIXME Remove following:
 786     delete textLength;     textLength = 0;
 787     delete textStartPos;   textStartPos = 0;
 788
 789     this->maxTextLength = 0;
 790     this->numberOfTexts = 0;
 791     this->numberOfAllTexts = 0;
 792     this->samplerate = samplerate;
 793     this->n = 0;
 794
 795     uchar verFlag = 0;
 796     if (std::fread(&verFlag, 1, 1, file) != 1)
 797         throw std::runtime_error("CSA::Load(): file read error (version flag).");
 798     if (verFlag != CSA::versionFlag)
 799         throw std::runtime_error("CSA::Load(): invalid save file version.");
 800
 801     if (std::fread(&(this->n), sizeof(TextPosition), 1, file) != 1)
 802         throw std::runtime_error("CSA::Load(): file read error (n).");
 803     if (std::fread(&samplerate, sizeof(unsigned), 1, file) != 1)
 804         throw std::runtime_error("CSA::Load(): file read error (samplerate).");
 805 // FIXME samplerate can not be changed during load.
 806 //    if (this->samplerate == 0)
 807 //        this->samplerate = samplerate;
 808
 809     for(ulong i = 0; i < 256; ++i)
 810         if (std::fread(this->C + i, sizeof(unsigned), 1, file) != 1)
 811             throw std::runtime_error("CSA::Load(): file read error (C table).");
 812
 813     if (std::fread(&(this->bwtEndPos), sizeof(TextPosition), 1, file) != 1)
 814         throw std::runtime_error("CSA::Load(): file read error (bwt end position).");
 815
 816     alphabetrank = static_sequence::load(file);
 817     sampled = new BSGAP(file);
 818     suffixes = new BlockArray(file);
 819     suffixDocId = new BlockArray(file);
 820
 821     if (std::fread(&(this->numberOfTexts), sizeof(unsigned), 1, file) != 1)
 822         throw std::runtime_error("CSA::Load(): file read error (numberOfTexts).");
 823     if (std::fread(&(this->numberOfAllTexts), sizeof(unsigned), 1, file) != 1)
 824         throw std::runtime_error("CSA::Load(): file read error (numberOfAllTexts).");
 825     if (std::fread(&(this->maxTextLength), sizeof(ulong), 1, file) != 1)
 826         throw std::runtime_error("CSA::Load(): file read error (maxTextLength).");
 827
 828     endmarkerDocId = new BlockArray(file);
 829     emptyTextRank = new BSGAP(file);
 830
 831     // FIXME Construct data structures with new samplerate
 832     //maketables();  // FIXME: this will redo text length tables
 833 }
 834
 835
 836
 837 /**
 838  * Rest of the functions follow...
 839  */
 840
 841 /*
 842  * Not supported
 843 uchar * CSA::Substring(TextPosition i, TextPosition l) const
 844 {
 845     uchar *result = new uchar[l + 1];
 846     if (l == 0)
 847     {
 848         result[0] = 0u;
 849         return result;
 850     }
 851
 852     TextPosition dist;
 853     TextPosition k = i + l - 1;
 854     // Check for end of the string
 855     if (k > n - 1)
 856     {
 857         l -= k - n + 1;
 858         k = n - 1;
 859     }
 860
 861     TextPosition skip = samplerate - k % samplerate - 1;
 862     TextPosition j;
 863     if (k / samplerate + 1 >= n / samplerate)
 864     {
 865         j = bwtEndPos;
 866         skip = n - k - 1;
 867     }
 868     else
 869     {
 870         j = (*positions)[k/samplerate+1];
 871         //cout << samplerate << ' ' << j << '\n';
 872     }
 873
 874     for (dist = 0; dist < skip + l; dist++)
 875     {
 876         int c = alphabetrank->access(j); //charAtPos(j, &alphabetrank_tmp);
 877         if (c == '\0')
 878         {
 879             // Rank among the end-markers in BWT
 880             unsigned endmarkerRank = alphabetrank->rank(0, j) - 1;
 881             j = (*endmarkerDocId)[endmarkerRank]; // LF-mapping for end-marker
 882         }
 883         else
 884             j = C[c]+alphabetrank->rank(c,j)-1; // LF-mapping
 885         if (dist >= skip)
 886             result[l + skip - dist - 1] = c;
 887     }
 888     result[l] = 0u;
 889     return result;
 890 }
 891 */
 892   /*TextPosition CSA::inverse(TextPosition i)
 893 {
 894     TextPosition skip = samplerate - i % samplerate;
 895     TextPosition j;
 896     if (i / samplerate + 1 >= n / samplerate)
 897     {
 898         j = bwtEndPos;
 899         skip = n - i;
 900     }
 901     else
 902     {
 903         j = (*positions)[i/samplerate+1];
 904         //cout << samplerate << ' ' << j << '\n';
 905     }
 906
 907     while (skip > 0)
 908     {
 909         int c = alphabetrank->charAtPos(j);
 910         j = C[c]+alphabetrank->rank(c,j)-1; // LF-mapping
 911         skip --;
 912     }
 913     return j;
 914     }*/
 915
 916 ulong CSA::Search(uchar const * pattern, TextPosition m, TextPosition *spResult, TextPosition *epResult) const
 917 {
 918     // use the FM-search replacing function Occ(c,1,i) with alphabetrank->rank(c,i)
 919     int c = (int)pattern[m-1];
 920     TextPosition i=m-1;
 921     TextPosition sp = C[c];
 922     TextPosition ep = C[c+1]-1;
 923     while (sp<=ep && i>=1)
 924     {
 925 //         printf("i = %lu, c = %c, sp = %lu, ep = %lu\n", i, pattern[i], sp, ep);
 926         c = (int)pattern[--i];
 927         sp = C[c]+alphabetrank->rank(c,sp-1);
 928         ep = C[c]+alphabetrank->rank(c,ep)-1;
 929     }
 930     *spResult = sp;
 931     *epResult = ep;
 932     if (sp<=ep)
 933         return ep - sp + 1;
 934     else
 935         return 0;
 936 }
 937
 938    /*TextPosition CSA::LF(uchar c, TextPosition &sp, TextPosition &ep)
 939 {
 940     sp = C[(int)c]+alphabetrank->rank(c,sp-1);
 941     ep = C[(int)c]+alphabetrank->rank(c,ep)-1;
 942
 943     if (sp<=ep)
 944         return ep - sp + 1;
 945     else
 946         return 0;
 947         }*/
 948
 949 CSA::TextPosition CSA::Lookup(TextPosition i) const
 950 {
 951     TextPosition dist=0;
 952     while (!sampled->IsBitSet(i))
 953     {
 954
 955         int c = alphabetrank->access(i);
 956         if (c == '\0')
 957         {
 958             // End-markers are sampled
 959             unsigned endmarkerRank = alphabetrank->rank(0, i) - 1;
 960             DocId docId = (*endmarkerDocId)[endmarkerRank];
 961             return (*textStartPos)[docId] + dist;
 962         }
 963         i = C[c]+alphabetrank->rank(c,i)-1;
 964         ++dist;
 965     }
 966
 967     return (*suffixes)[sampled->rank(i)-1]+dist;
 968 }
 969
 970 CSA::~CSA() {
 971     delete dynFMI;
 972     delete alphabetrank;
 973     delete sampled;
 974     delete suffixes;
 975     delete suffixDocId;
 976     delete positions;
 977     delete [] codetable;
 978
 979     delete endmarkerDocId;
 980     delete textLength;
 981     delete textStartPos;
 982     delete emptyTextRank;
 983 }
 984
 985 void CSA::makewavelet(uchar *bwt)
 986 {
 987     ulong i, min = 0,
 988              max;
 989     for (i=0;i<256;i++)
 990         C[i]=0;
 991     for (i=0;i<n;++i)
 992         C[(int)bwt[i]]++;
 993     for (i=0;i<256;i++)
 994         if (C[i]>0) {min = i; break;}
 995     for (i=255;i>=min;--i)
 996         if (C[i]>0) {max = i; break;}
 997
 998     // Print frequencies
 999     /*for(i = 0; i < 256; i++)
1000         if (C[i]>0) printf("C[%lu] = %lu\n", i, C[i]);
1001         fflush(stdout);*/
1002
1003     ulong prev=C[0], temp;
1004     C[0]=0;
1005     for (i=1;i<256;i++) {
1006         temp = C[i];
1007         C[i]=C[i-1]+prev;
1008         prev = temp;
1009     }
1010 //    this->codetable = node::makecodetable(bwt,n);
1011 //    alphabetrank = new THuffAlphabetRank(bwt,n, this->codetable,0);
1012 //    delete [] bwt;
1013     //alphabetrank = new RLWaveletTree(bwt, n); // Deletes bwt!
1014
1015     uint *text = new uint[n];
1016     for (i = 0; i < n; ++i) // Silly
1017         text[i] = bwt[i];
1018     delete [] bwt;
1019     bwt = 0;
1020
1021     alphabet_mapper * am = new alphabet_mapper_none();
1022     static_bitsequence_builder * bmb = new static_bitsequence_builder_brw32(16); // FIXME samplerate?
1023
1024     //cout << "Building Huffman table..."; cout.flush();
1025
1026     wt_coder * wtc = new wt_coder_huff(text,n,am);
1027
1028     //cout << "done" << endl; cout.flush();
1029     //cout << "Building static_sequence..."; cout.flush();
1030
1031     alphabetrank = new static_sequence_wvtree(text,n,wtc,bmb,am);
1032     delete [] text;
1033     text = 0;
1034
1035 /*    for (i = 0; i < n; ++i)
1036     {
1037         uchar c = alphabetrank->charAtPos(i);
1038         TextPosition j = C[c]+alphabetrank->rank(c, i)-1;
1039         printf("LF[%lu] = %lu\n", i, j);
1040         }*/
1041 }
1042
1043 void CSA::maketables()
1044 {
1045     ulong sampleLength = (n%samplerate==0) ? n/samplerate : n/samplerate+1;
1046
1047     ulong *sampledpositions = new ulong[n/W+1];
1048     unsigned ceilLog2n = Tools::CeilLog2(n);
1049     positions = new BlockArray(sampleLength, ceilLog2n);
1050     BlockArray* tmpSuffix = new BlockArray(sampleLength, ceilLog2n);
1051
1052     // Mapping from end-markers to doc ID's:
1053     endmarkerDocId = new BlockArray(numberOfTexts, Tools::CeilLog2(numberOfTexts));
1054
1055     ulong i,j=0;
1056     for (i=0;i<n/W+1;i++)
1057         sampledpositions[i]=0lu;
1058
1059     ulong x,p=bwtEndPos;
1060     ulong sampleCount = 0;
1061     // Keeping track of text position of end-markers seen
1062     ulong posOfSuccEndmarker = n;
1063     DocId textId = numberOfTexts;
1064     ulong ulongmax = 0;
1065     ulongmax--;
1066
1067     //positions:
1068     for (i=n-1;i<ulongmax;i--) { // TODO bad solution with ulongmax?
1069       // i substitutes SA->GetPos(i)
1070         x=(i==n-1)?0:i+1;
1071
1072         if (x % samplerate == 0 && posOfSuccEndmarker - x > samplerate) {
1073             Tools::SetField(sampledpositions,1,p,1);
1074             (*positions)[sampleCount] = p;
1075             (*tmpSuffix)[sampleCount] = x; // FIXME remove
1076             sampleCount ++;
1077         }
1078
1079         uchar c = alphabetrank->access(p);
1080         if (c == '\0')
1081         {
1082             --textId;
1083
1084             // Record the order of end-markers in BWT:
1085             ulong endmarkerRank = alphabetrank->rank(0, p) - 1;
1086             (*endmarkerDocId)[endmarkerRank] = textId;
1087
1088             // Store text length and text start position:
1089             if (textId < (DocId)numberOfTexts - 1)
1090             {
1091                 (*textLength)[textId + 1] = posOfSuccEndmarker - x;
1092                 (*textStartPos)[textId + 1] = x;  // x is the position of end-marker.
1093                 posOfSuccEndmarker = x;
1094             }
1095
1096             // LF-mapping from '\0' does not work with this (pseudo) BWT (see details from Wolfgang's thesis).
1097             p = textId; // Correct LF-mapping to the last char of the previous text.
1098         }
1099         else
1100             p = C[c]+alphabetrank->rank(c, p)-1;
1101     }
1102     assert(textId == 0);
1103
1104     /*i = 0;
1105     for (map<ulong, pair<DocId, ulong> >::iterator it = endmarkers.begin(); it != endmarkers.end(); ++it, ++i)
1106     {
1107         int docc = (*endmarkerDocId)[i];
1108         ulong poss = (*endmarkerPos)[i];
1109         printf("endm[%u] = %lu (text pos: %lu) (recorded: %d, %lu)\n", (it->second).first, it->first, (it->second).second, docc, poss);
1110         }*/
1111 /*
1112     for (i = 0; i < numberOfTexts; ++ i)
1113     {
1114         //std::cout << "textlength = " << dynTextLength[i].first << " vs " << (*textLength)[i] << ", textStartPos = " << dynTextLength[i].second << " vs " << (*textStartPos)[i] << std::endl;
1115         assert(dynTextLength[i].first == (*textLength)[i]);
1116         assert(dynTextLength[i].second == (*textStartPos)[i]);
1117     }
1118 */
1119     sampled = new BSGAP(sampledpositions,n,true);
1120     sampleLength = sampled->rank(n-1);
1121     assert(sampleCount == sampleLength);
1122 //    std::cout << ";sampleLength;" << sampleLength << std::endl;
1123     // Suffixes == offset from text start position
1124     suffixes = new BlockArray(sampleLength, Tools::CeilLog2(maxTextLength));
1125     suffixDocId = new BlockArray(sampleLength, Tools::CeilLog2(numberOfTexts));
1126
1127    //suffixes:
1128     for(i=0; i<sampleLength; i++) {
1129         assert((*positions)[i] < n);
1130         j = sampled->rank((*positions)[i]);
1131         if (j==0) j=sampleLength;
1132         TextPosition textPos = (*tmpSuffix)[i]; //(i*samplerate==n)?0:i*samplerate;
1133         (*suffixDocId)[j-1] = DocIdAtTextPos(textPos); // (*suffixes)[j-1]);
1134
1135         assert((unsigned)DocIdAtTextPos(textPos) < numberOfTexts);
1136         assert((*suffixDocId)[j-1] < numberOfTexts);
1137         // calculate offset from text start:
1138         (*suffixes)[j-1] = textPos - (*textStartPos)[(*suffixDocId)[j-1]];
1139     }
1140     // FIXME Temp, remove
1141     delete tmpSuffix;
1142     tmpSuffix=0;
1143     delete positions;
1144     positions =0;
1145     delete textLength;
1146     textLength = 0;
1147     delete textStartPos;
1148     textStartPos = 0;
1149 }
1150
1151
1152 /**
1153  * Finds document identifier for given text position
1154  *
1155  * Starting text position of the document is stored into second parameter.
1156  * Binary searching on text starting positions.
1157  */
1158 TextCollection::DocId CSA::DocIdAtTextPos(TextPosition i) const
1159 {
1160     assert(i < n);
1161
1162     DocId a = 0;
1163     DocId b = numberOfTexts - 1;
1164     while (a < b)
1165     {
1166         DocId c = a + (b - a)/2;
1167         if ((*textStartPos)[c] > i)
1168             b = c - 1;
1169         else if ((*textStartPos)[c+1] > i)
1170             return c;
1171         else
1172             a = c + 1;
1173     }
1174
1175     assert(a < (DocId)numberOfTexts);
1176     assert(i >= (*textStartPos)[a]);
1177     assert(i < (a == (DocId)numberOfTexts - 1 ? n : (*textStartPos)[a+1]));
1178     return a;
1179 }
1180
1181 CSA::TCodeEntry * CSA::node::makecodetable(uchar *text, TextPosition n)
1182 {
1183     TCodeEntry *result = new TCodeEntry[ 256 ];
1184
1185     count_chars( text, n, result );
1186     std::priority_queue< node, std::vector< node >, std::greater<node> > q;
1187 //
1188 // First I push all the leaf nodes into the queue
1189 //
1190     for ( unsigned int i = 0 ; i < 256 ; i++ )
1191         if ( result[ i ].count )
1192             q.push(node( i, result[ i ].count ) );
1193 //
1194 // This loop removes the two smallest nodes from the
1195 // queue.  It creates a new internal node that has
1196 // those two nodes as children. The new internal node
1197 // is then inserted into the priority queue.  When there
1198 // is only one node in the priority queue, the tree
1199 // is complete.
1200 //
1201
1202     while ( q.size() > 1 ) {
1203         node *child0 = new node( q.top() );
1204         q.pop();
1205         node *child1 = new node( q.top() );
1206         q.pop();
1207         q.push( node( child0, child1 ) );
1208     }
1209 //
1210 // Now I compute and return the codetable
1211 //
1212     q.top().maketable(0u,0u, result);
1213     q.pop();
1214     return result;
1215 }
1216
1217
1218 void CSA::node::maketable(unsigned code, unsigned bits, TCodeEntry *codetable) const
1219 {
1220     if ( child0 )
1221     {
1222         child0->maketable( SetBit(code,bits,0), bits+1, codetable );
1223         child1->maketable( SetBit(code,bits,1), bits+1, codetable );
1224         delete child0;
1225         delete child1;
1226     }
1227     else
1228     {
1229         codetable[value].code = code;
1230         codetable[value].bits = bits;
1231     }
1232 }
1233
1234 void CSA::node::count_chars(uchar *text, TextPosition n, TCodeEntry *counts )
1235 {
1236     ulong i;
1237     for (i = 0 ; i < 256 ; i++ )
1238         counts[ i ].count = 0;
1239     for (i=0; i<n; i++)
1240         counts[(int)text[i]].count++;
1241 }
1242
1243 unsigned CSA::node::SetBit(unsigned x, unsigned pos, unsigned bit) {
1244       return x | (bit << pos);
1245 }
1246