incbwt/rlcsa.cpp

   1 #include <algorithm>
   2 #include <cstdlib>
   3 #include <cstring>
   4 #include <fstream>
   5 #include <iostream>
   6 #include <list>
   7
   8 #include "rlcsa.h"
   9 #include "misc/utils.h"
  10 #include "qsufsort/qsufsort.h"
  11 #include "bits/vectors.h"
  12
  13
  14 #ifdef MULTITHREAD_SUPPORT
  15 #include <omp.h>
  16 #endif
  17
  18
  19
  20 namespace CSA
  21 {
  22
  23
  24 RLCSA::RLCSA(const std::string& base_name, bool print) :
  25   ok(false),
  26   sa_samples(0),
  27   end_points(0)
  28 {
  29   for(usint c = 0; c < CHARS; c++) { this->array[c] = 0; }
  30
  31   std::string array_name = base_name + ARRAY_EXTENSION;
  32   std::ifstream array_file(array_name.c_str(), std::ios_base::binary);
  33   if(!array_file)
  34   {
  35     std::cerr << "RLCSA: Error opening Psi array file!" << std::endl;
  36     return;
  37   }
  38
  39   usint distribution[CHARS];
  40   array_file.read((char*)distribution, CHARS * sizeof(usint));
  41   this->buildCharIndexes(distribution);
  42
  43 //  Parameters parameters;
  44 //  parameters.read(base_name + PARAMETERS_EXTENSION);
  45   for(usint c = 0; c < CHARS; c++)
  46   {
  47     if(!isEmpty(this->index_ranges[c])) { this->array[c] = new PsiVector(array_file); }
  48   }
  49
  50   this->end_points = new DeltaVector(array_file);
  51   this->number_of_sequences = this->end_points->getNumberOfItems();
  52
  53   array_file.read((char*)&(this->sample_rate), sizeof(this->sample_rate));
  54   array_file.close();
  55
  56   this->support_locate = true;// parameters.get(SUPPORT_LOCATE);
  57   this->support_display = true; //parameters.get(SUPPORT_DISPLAY);
  58
  59   if(this->support_locate || this->support_display)
  60   {
  61     std::string sa_sample_name = base_name + SA_SAMPLES_EXTENSION;
  62     std::ifstream sa_sample_file(sa_sample_name.c_str(), std::ios_base::binary);
  63     if(!sa_sample_file)
  64     {
  65       std::cerr << "RLCSA: Error opening suffix array sample file!" << std::endl;
  66       return;
  67     }
  68     this->sa_samples = new SASamples(sa_sample_file, this->sample_rate);
  69     sa_sample_file.close();
  70   }
  71
  72 //  if(print) { parameters.print(); }
  73
  74   this->ok = true;
  75 }
  76
  77 RLCSA::RLCSA(uchar* data, usint bytes, usint block_size, usint sa_sample_rate, bool multiple_sequences, bool delete_data) :
  78   ok(false),
  79   sa_samples(0),
  80   sample_rate(sa_sample_rate),
  81   end_points(0)
  82 {
  83   for(usint c = 0; c < CHARS; c++) { this->array[c] = 0; }
  84
  85   if(!data || bytes == 0)
  86   {
  87     std::cerr << "RLCSA: No input data given!" << std::endl;
  88     if(delete_data) { delete[] data; }
  89     return;
  90   }
  91   if(block_size < 2 * sizeof(usint) || block_size % sizeof(usint) != 0)
  92   {
  93     std::cerr << "RLCSA: Block size must be a multiple of " << sizeof(usint) << " bytes!" << std::endl;
  94     if(delete_data) { delete[] data; }
  95     return;
  96   }
  97
  98
  99   // Do we store SA samples?
 100   if(this->sample_rate > 0)
 101   {
 102     this->support_locate = this->support_display = true;
 103   }
 104   else
 105   {
 106     this->support_locate = this->support_display = false;
 107     this->sample_rate = 1;
 108   }
 109
 110
 111   // Determine the number of sequences and mark their end points.
 112   if(multiple_sequences)
 113   {
 114     DeltaEncoder endings(RLCSA::ENDPOINT_BLOCK_SIZE);
 115     this->number_of_sequences = 0;
 116     usint marker = 0;
 117     usint padding = 0, chars_encountered = 0;
 118
 119     for(usint i = 0; i < bytes; i++)
 120     {
 121       if(data[i] == 0)
 122       {
 123         if(i == marker) { break; }  // Empty sequence.
 124         this->number_of_sequences++;
 125         marker = i + 1;
 126
 127         usint pos = chars_encountered + padding - 1;
 128         endings.setBit(pos);
 129         padding = ((pos + this->sample_rate) / this->sample_rate) * this->sample_rate - chars_encountered;
 130       }
 131       else
 132       {
 133         chars_encountered++;
 134       }
 135     }
 136
 137     if(this->number_of_sequences == 0 || marker != bytes)
 138     {
 139       std::cerr << "RLCSA: Collection must consist of 0-terminated nonempty sequences!" << std::endl;
 140       if(delete_data) { delete[] data; }
 141       return;
 142     }
 143     this->end_points = new DeltaVector(endings, chars_encountered + padding);
 144   }
 145   else
 146   {
 147     this->number_of_sequences = 1;
 148     DeltaEncoder endings(RLCSA::ENDPOINT_BLOCK_SIZE, RLCSA::ENDPOINT_BLOCK_SIZE);
 149     usint pos = bytes - 1;
 150     endings.setBit(pos);
 151     usint padding = ((pos + this->sample_rate) / this->sample_rate) * this->sample_rate - bytes;
 152     this->end_points = new DeltaVector(endings, bytes + padding);
 153   }
 154
 155
 156   // Build character tables etc.
 157   usint distribution[CHARS];
 158   sint low = CHARS, high = 0;
 159   for(usint c = 0; c < CHARS; c++) { distribution[c] = 0; }
 160   for(usint i = 0; i < bytes; i++)
 161   {
 162     if(data[i] < low) { low = data[i]; }
 163     if(data[i] > high) { high = data[i]; }
 164     distribution[(usint)data[i]]++;
 165   }
 166   if(multiple_sequences)
 167   {
 168     distribution[0] = 0;
 169     low = 0;
 170     high = high + this->number_of_sequences;
 171   }
 172   this->buildCharIndexes(distribution);
 173
 174
 175   // Build suffix array.
 176   int* inverse = new int[bytes + 1];
 177   if(multiple_sequences)
 178   {
 179     int zeros = 0;
 180     for(usint i = 0; i < bytes; i++)
 181     {
 182       if(data[i] == 0)
 183       {
 184         inverse[i] = zeros;
 185         zeros++;
 186       }
 187       else
 188       {
 189         inverse[i] = (int)(data[i] + this->number_of_sequences);
 190       }
 191     }
 192   }
 193   else
 194   {
 195     for(usint i = 0; i < bytes; i++) { inverse[i] = (int)data[i]; }
 196   }
 197   if(delete_data) { delete[] data; }
 198   int* sa = new int[bytes + 1];
 199   suffixsort(inverse, sa, bytes, high + 1, low);
 200
 201
 202   // Sample SA.
 203   if(this->support_locate || this->support_display)
 204   {
 205     if(multiple_sequences)
 206     {
 207       this->sa_samples = new SASamples((uint*)inverse, this->end_points, this->data_size, this->sample_rate);
 208     }
 209     else
 210     {
 211       this->sa_samples = new SASamples((uint*)(sa + 1), this->data_size, this->sample_rate);
 212     }
 213   }
 214
 215
 216   // Build Psi.
 217   for(usint i = 0; i <= bytes; i++)
 218   {
 219     sa[i] = inverse[(sa[i] + 1) % (bytes + 1)];
 220   }
 221   delete[] inverse;
 222
 223
 224   // Build RLCSA.
 225   usint incr = (multiple_sequences ? this->number_of_sequences + 1 : 1);  // No e_of_s markers in SA.
 226   usint decr = (multiple_sequences ? 1 : 0);  // Ignore the last e_of_s marker if multiple sequences.
 227   for(usint c = 0; c < CHARS; c++)
 228   {
 229     if(distribution[c] == 0) { this->array[c] = 0; continue; }
 230
 231     uint* curr = (uint*)(sa + index_ranges[c].first + incr);
 232     uint* limit = (uint*)(sa + index_ranges[c].second + incr + 1);
 233     PsiVector::Encoder encoder(block_size);
 234     pair_type run(*curr, 1);
 235     curr++;
 236
 237     for(; curr < limit; curr++)
 238     {
 239       if(*curr == run.first + run.second) { run.second++; }
 240       else
 241       {
 242         encoder.setRun(run.first - decr, run.second);
 243         run = pair_type(*curr, 1);
 244       }
 245     }
 246     encoder.setRun(run.first - decr, run.second);
 247
 248     this->array[c] = new PsiVector(encoder, this->data_size + incr - decr);
 249   }
 250   delete[] sa;
 251
 252
 253   this->ok = true;
 254 }
 255
 256 RLCSA::RLCSA(RLCSA& index, RLCSA& increment, usint* positions, usint block_size, usint threads) :
 257   ok(false),
 258   sa_samples(0),
 259   end_points(0)
 260 {
 261   for(usint c = 0; c < CHARS; c++) { this->array[c] = 0; }
 262
 263   if(!index.isOk() || !increment.isOk())
 264   {
 265     return; // Fail silently. Actual error has already been reported.
 266   }
 267   if(positions == 0)
 268   {
 269     std::cerr << "RLCSA: Positions for insertions not available!" << std::endl;
 270     return;
 271   }
 272   if(index.sample_rate != increment.sample_rate)
 273   {
 274     std::cerr << "RLCSA: Cannot combine indexes with different sample rates!" << std::endl;
 275     return;
 276   }
 277
 278   if(index.sa_samples == 0 || increment.sa_samples == 0)
 279   {
 280     this->support_locate = this->support_display = false;
 281   }
 282   else
 283   {
 284     this->support_locate = this->support_display = true;
 285   }
 286
 287
 288   // Build character tables etc.
 289   usint distribution[CHARS];
 290   for(usint c = 0; c < CHARS; c++)
 291   {
 292     distribution[c] = length(index.index_ranges[c]) + length(increment.index_ranges[c]);
 293   }
 294   this->buildCharIndexes(distribution);
 295   this->sample_rate = index.sample_rate;
 296   this->number_of_sequences = index.number_of_sequences + increment.number_of_sequences;
 297
 298   // Merge end points, SA samples, and Psi.
 299   usint psi_size = this->data_size + this->number_of_sequences;
 300   bool should_be_ok = true;
 301
 302   #ifdef MULTITHREAD_SUPPORT
 303   omp_set_num_threads(threads);
 304   #pragma omp parallel
 305   {
 306     #pragma omp for schedule(dynamic, 1)
 307   #endif
 308     for(int c = -2; c < (int)CHARS; c++)
 309     {
 310       if(c == -2)      { this->mergeEndPoints(index, increment); }
 311       else if(c == -1) { this->mergeSamples(index, increment, positions);  }
 312       else if(distribution[c] != 0)
 313       {
 314         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);
 315         index.array[c] = 0;
 316         increment.array[c] = 0;
 317
 318         if(this->array[c] == 0)
 319         {
 320           std::cerr << "RLCSA: Merge failed for vectors " << c << "!" << std::endl;
 321           should_be_ok = false;
 322         }
 323       }
 324     }
 325   #ifdef MULTITHREAD_SUPPORT
 326   }
 327   #endif
 328
 329   this->ok = should_be_ok;
 330 }
 331
 332 RLCSA::~RLCSA()
 333 {
 334   for(usint c = 0; c < CHARS; c++) { delete this->array[c]; }
 335   delete this->sa_samples;
 336   delete this->end_points;
 337 }
 338
 339 //--------------------------------------------------------------------------
 340
 341 void
 342 RLCSA::writeTo(const std::string& base_name) const
 343 {
 344   std::string array_name = base_name + ARRAY_EXTENSION;
 345   std::ofstream array_file(array_name.c_str(), std::ios_base::binary);
 346   if(!array_file)
 347   {
 348     std::cerr << "RLCSA: Error creating Psi array file!" << std::endl;
 349     return;
 350   }
 351
 352   usint distribution[CHARS];
 353   for(usint c = 0; c < CHARS; c++)
 354   {
 355     distribution[c] = length(this->index_ranges[c]);
 356   }
 357   array_file.write((char*)distribution, CHARS * sizeof(usint));
 358   for(usint c = 0; c < CHARS; c++)
 359   {
 360     if(this->array[c] != 0)
 361     {
 362       this->array[c]->writeTo(array_file);
 363     }
 364   }
 365
 366   this->end_points->writeTo(array_file);
 367   array_file.write((char*)&(this->sample_rate), sizeof(this->sample_rate));
 368   array_file.close();
 369
 370   if(this->support_locate || this->support_display)
 371   {
 372     std::string sa_sample_name = base_name + SA_SAMPLES_EXTENSION;
 373     std::ofstream sa_sample_file(sa_sample_name.c_str(), std::ios_base::binary);
 374     if(!sa_sample_file)
 375     {
 376       std::cerr << "RLCSA: Error creating suffix array sample file!" << std::endl;
 377       return;
 378     }
 379
 380     this->sa_samples->writeTo(sa_sample_file);
 381     sa_sample_file.close();
 382   }
 383 }
 384
 385 bool
 386 RLCSA::isOk() const
 387 {
 388   return this->ok;
 389 }
 390
 391 //--------------------------------------------------------------------------
 392
 393 pair_type
 394 RLCSA::count(const std::string& pattern) const
 395 {
 396   if(pattern.length() == 0) { return pair_type(0, this->data_size - 1); }
 397
 398   std::string::const_reverse_iterator iter = pattern.rbegin();
 399
 400   pair_type index_range = this->index_ranges[(uchar)*iter];
 401   if(isEmpty(index_range)) { return EMPTY_PAIR; }
 402   index_range.first += this->number_of_sequences;
 403   index_range.second += this->number_of_sequences;
 404
 405   for(++iter; iter != pattern.rend(); ++iter)
 406   {
 407     index_range = this->backwardSearchStep(index_range, (uchar)*iter);
 408     if(isEmpty(index_range)) { return EMPTY_PAIR; }
 409   }
 410
 411   // Suffix array indexes are 0-based.
 412   index_range.first -= this->number_of_sequences;
 413   index_range.second -= this->number_of_sequences;
 414
 415   return index_range;
 416 }
 417
 418 //--------------------------------------------------------------------------
 419
 420 void
 421 RLCSA::reportPositions(uchar* data, usint length, usint* positions) const
 422 {
 423   if(data == 0 || length == 0 || positions == 0) { return; }
 424
 425   PsiVector::Iterator** iters = this->getIterators();
 426
 427   usint current = this->number_of_sequences - 1;
 428   positions[length] = current; // "immediately after current"
 429   for(sint i = (sint)(length - 1); i >= 0; i--)
 430   {
 431     usint c = (usint)data[i];
 432     if(this->array[c] != 0)
 433     {
 434       current = this->LF(current, c, *(iters[c]));
 435     }
 436     else
 437     {
 438       if(c < this->text_chars[0]) // No previous characters either.
 439       {
 440         current = this->number_of_sequences - 1;
 441       }
 442       else
 443       {
 444         current = this->index_ranges[c].first - 1 + this->number_of_sequences;
 445       }
 446     }
 447     positions[i] = current; // "immediately after current"
 448   }
 449
 450   this->deleteIterators(iters);
 451 }
 452
 453 //--------------------------------------------------------------------------
 454
 455 usint*
 456 RLCSA::locate(pair_type range) const
 457 {
 458   if(!(this->support_locate) || isEmpty(range) || range.second >= this->data_size) { return 0; }
 459
 460   usint* data = new usint[length(range)];
 461   this->locateUnsafe(range, data);
 462
 463   return data;
 464 }
 465
 466 usint*
 467 RLCSA::locate(pair_type range, usint* data) const
 468 {
 469   if(!(this->support_locate) || isEmpty(range) || range.second >= this->data_size || data == 0) { return 0; }
 470   this->locateUnsafe(range, data);
 471   return data;
 472 }
 473
 474 void
 475 RLCSA::locateUnsafe(pair_type range, usint* data) const
 476 {
 477   usint items = length(range);
 478   usint* offsets = new usint[items];
 479   bool* finished = new bool[items];  // FIXME This could be more space efficient...
 480
 481   PsiVector::Iterator** iters = this->getIterators();
 482
 483   for(usint i = 0, j = range.first; i < items; i++, j++)
 484   {
 485     data[i] = j + this->number_of_sequences;
 486     offsets[i] = 0;
 487     finished[i] = false;
 488   }
 489
 490   bool found = false;
 491   while(!found)
 492   {
 493     found = true;
 494     pair_type run = EMPTY_PAIR;
 495     for(usint i = 0; i < items; i++)
 496     {
 497       if(finished[i])
 498       {
 499         continue; // The run might continue after this.
 500       }
 501       else if(isEmpty(run))
 502       {
 503         run = pair_type(i, i);
 504       }
 505       else if(data[i] - data[run.first] == i - run.first)
 506       {
 507         run.second = i;
 508       }
 509       else
 510       {
 511         found &= this->processRun(run, data, offsets, finished, iters);
 512         run = pair_type(i, i);
 513       }
 514     }
 515     if(!isEmpty(run)) { found &= this->processRun(run, data, offsets, finished, iters); }
 516   }
 517
 518   this->deleteIterators(iters);
 519   delete[] offsets;
 520   delete[] finished;
 521 }
 522
 523 bool
 524 RLCSA::processRun(pair_type run, usint* data, usint* offsets, bool* finished, PsiVector::Iterator** iters) const
 525 {
 526   bool found = true;
 527   usint run_start = 0, run_left = 0;
 528   pair_type next_sample = pair_type(0, 0);
 529
 530   for(usint i = run.first; i <= run.second; i++)
 531   {
 532     if(finished[i])
 533     {
 534       if(run_left > 0) { run_left--; }
 535       continue;
 536     }
 537     if(data[i] < this->number_of_sequences) // Implicit sample here.
 538     {
 539       DeltaVector::Iterator iter(*(this->end_points));
 540       data[i] = iter.select(data[i]) + 1 - offsets[i];
 541       finished[i] = true;
 542       if(run_left > 0) { run_left--; }
 543       continue;
 544     }
 545     if(next_sample.first < data[i]) // Need another sample.
 546     {
 547       next_sample = this->sa_samples->getFirstSampleAfter(data[i] - this->number_of_sequences);
 548       next_sample.first += this->number_of_sequences;
 549     }
 550     if(data[i] < next_sample.first) // No sample found for current position.
 551     {
 552       if(run_left > 0)
 553       {
 554         data[i] = data[run_start] + i - run_start;
 555         run_left--;
 556       }
 557       else
 558       {
 559         pair_type value = this->psi(data[i] - this->number_of_sequences, run.second - i, iters);
 560         data[i] = value.first;
 561         run_left = value.second;
 562         run_start = i;
 563       }
 564       offsets[i]++;
 565       found = false;
 566     }
 567     else  // Sampled position found.
 568     {
 569       data[i] = this->sa_samples->getSample(next_sample.second) - offsets[i];
 570       finished[i] = true;
 571       if(run_left > 0) { run_left--; }
 572     }
 573   }
 574   return found;
 575 }
 576
 577 usint
 578 RLCSA::locate(usint index) const
 579 {
 580   if(!(this->support_locate) || index >= this->data_size) { return this->data_size; }
 581
 582   usint offset = 0;
 583   index += this->number_of_sequences;
 584   while(true)
 585   {
 586     if(index < this->number_of_sequences) // Implicit sample here
 587     {
 588       DeltaVector::Iterator iter(*(this->end_points));
 589       return iter.select(index) + 1 - offset;
 590     }
 591     pair_type next_sample = this->sa_samples->getFirstSampleAfter(index - this->number_of_sequences);
 592     next_sample.first += this->number_of_sequences;
 593     if(next_sample.first == index)
 594     {
 595       return this->sa_samples->getSample(next_sample.second) - offset;
 596     }
 597     index = this->psi(index - this->number_of_sequences);
 598     offset++;
 599   }
 600 }
 601
 602 //--------------------------------------------------------------------------
 603
 604 uchar*
 605 RLCSA::display(usint sequence) const
 606 {
 607   if(!(this->support_display)) { return 0; }
 608
 609   pair_type seq_range = this->getSequenceRange(sequence);
 610   if(isEmpty(seq_range)) { return 0; }
 611
 612   uchar* data = new uchar[length(seq_range)+1];
 613   this->displayUnsafe(seq_range, data);
 614   data[length(seq_range)] = 0;
 615
 616   return data;
 617 }
 618
 619 uchar*
 620 RLCSA::display(usint sequence, pair_type range) const
 621 {
 622   if(!(this->support_display) || isEmpty(range)) { return 0; }
 623
 624   pair_type seq_range = this->getSequenceRange(sequence);
 625   if(isEmpty(seq_range)) { return 0; }
 626
 627   range.first += seq_range.first; range.second += seq_range.first;
 628   if(range.second > seq_range.second) { return 0; }
 629
 630   uchar* data = new uchar[length(range)];
 631   this->displayUnsafe(range, data);
 632
 633   return data;
 634 }
 635
 636 uchar*
 637 RLCSA::display(usint sequence, pair_type range, uchar* data) const
 638 {
 639   if(!(this->support_display) || isEmpty(range) || data == 0) { return 0; }
 640
 641   pair_type seq_range = this->getSequenceRange(sequence);
 642   if(isEmpty(seq_range)) { return 0; }
 643
 644   range.first += seq_range.first; range.second += seq_range.first;
 645   if(range.second > seq_range.second) { return 0; }
 646
 647   this->displayUnsafe(range, data);
 648   return data;
 649 }
 650
 651 uchar*
 652 RLCSA::display(usint position, usint len, usint context, usint& result_length) const
 653 {
 654   if(!(this->support_display)) { return 0; }
 655
 656   pair_type range = this->getSequenceRangeForPosition(position);
 657   if(isEmpty(range)) { return 0; }
 658
 659   range.first = position - std::min(context, position - range.first);
 660   range.second = std::min(range.second, position + len + context - 1);
 661   result_length = length(range);
 662   if(isEmpty(range)) { return 0; }
 663
 664   uchar* data = new uchar[length(range)];
 665   this->displayUnsafe(range, data);
 666
 667   return data;
 668 }
 669
 670 usint
 671 RLCSA::displayPrefix(usint sequence, usint len, uchar* data) const
 672 {
 673   if(!(this->support_display) || len == 0 || data == 0) { return 0; }
 674
 675   pair_type seq_range = this->getSequenceRange(sequence);
 676   if(isEmpty(seq_range)) { return 0; }
 677
 678   pair_type range(seq_range.first, std::min(seq_range.second, seq_range.first + len - 1));
 679
 680   this->displayUnsafe(range, data);
 681   return length(range);
 682 }
 683
 684 usint
 685 RLCSA::displayFromPosition(usint index, usint max_len, uchar* data) const
 686 {
 687   if(max_len == 0 || data == 0 || index >= this->data_size) { return 0; }
 688
 689   for(usint i = 0; i < max_len; i++)
 690   {
 691     data[i] = this->getCharacter(index);
 692     index = this->psiUnsafe(index, data[i]);
 693     if(index < this->number_of_sequences) { return i + 1; }
 694     index -= this->number_of_sequences;
 695   }
 696
 697   return max_len;
 698 }
 699
 700 void
 701 RLCSA::displayUnsafe(pair_type range, uchar* data) const
 702 {
 703   usint i = range.first - range.first % this->sa_samples->getSampleRate();
 704
 705   usint pos = this->sa_samples->inverseSA(i);
 706
 707   if(length(range) >= 1024)
 708   {
 709     PsiVector::Iterator** iters = this->getIterators();
 710     for(; i < range.first; i++)
 711     {
 712       pos = this->psi(pos, iters) - this->number_of_sequences;
 713     }
 714     for(; i <= range.second; i++)
 715     {
 716       usint c = this->getCharacter(pos);
 717       data[i - range.first] = c;
 718       pos = this->psiUnsafe(pos, c, *(iters[c])) - this->number_of_sequences;
 719     }
 720     this->deleteIterators(iters);
 721   }
 722   else
 723   {
 724     for(; i < range.first; i++)
 725     {
 726       pos = this->psi(pos) - this->number_of_sequences;
 727     }
 728     for(; i <= range.second; i++)
 729     {
 730       usint c = this->getCharacter(pos);
 731       data[i - range.first] = c;
 732       pos = this->psiUnsafe(pos, c) - this->number_of_sequences;
 733     }
 734   }
 735 }
 736
 737 //--------------------------------------------------------------------------
 738
 739 pair_type
 740 RLCSA::getSequenceRange(usint number) const
 741 {
 742   if(number >= this->number_of_sequences) { return EMPTY_PAIR; }
 743
 744   pair_type result;
 745   DeltaVector::Iterator iter(*(this->end_points));
 746   if(number == 0)
 747   {
 748     result.first = 0;
 749     result.second = iter.select(number);
 750   }
 751   else
 752   {
 753     result.first = nextMultipleOf(this->sample_rate, iter.select(number - 1));
 754     result.second = iter.selectNext();
 755   }
 756
 757   return result;
 758 }
 759
 760 usint
 761 RLCSA::getSequenceForPosition(usint value) const
 762 {
 763   if(value == 0) { return 0; }
 764   DeltaVector::Iterator iter(*(this->end_points));
 765   return iter.rank(value - 1);
 766 }
 767
 768 pair_type
 769 RLCSA::getSequenceRangeForPosition(usint value) const
 770 {
 771   return this->getSequenceRange(this->getSequenceForPosition(value));
 772 }
 773
 774 usint*
 775 RLCSA::getSequenceForPosition(usint* values, usint len) const
 776 {
 777   if(values == 0) { return 0; }
 778
 779   DeltaVector::Iterator iter(*(this->end_points));
 780   for(usint i = 0; i < len; i++)
 781   {
 782     if(values[i] > 0) { values[i] = iter.rank(values[i] - 1); }
 783   }
 784
 785   return values;
 786 }
 787
 788 pair_type
 789 RLCSA::getRelativePosition(usint value) const
 790 {
 791   DeltaVector::Iterator iter(*(this->end_points));
 792   pair_type result(0, value);
 793
 794   if(value > 0) { result.first = iter.rank(value - 1); }
 795   if(result.first > 0)
 796   {
 797     result.second -= nextMultipleOf(this->sample_rate, iter.select(result.first - 1));
 798   }
 799
 800   return result;
 801 }
 802
 803 //--------------------------------------------------------------------------
 804
 805 uchar*
 806 RLCSA::readBWT() const
 807 {
 808   return this->readBWT(pair_type(0, this->data_size + this->number_of_sequences - 1));
 809 }
 810
 811 uchar*
 812 RLCSA::readBWT(pair_type range) const
 813 {
 814   if(isEmpty(range) || range.second >= this->data_size + this->number_of_sequences) { return 0; }
 815
 816   usint n = length(range);
 817
 818   uchar* bwt = new uchar[n];
 819   memset(bwt, 0, n);
 820
 821   for(usint c = 0; c < CHARS; c++)
 822   {
 823     if(this->array[c] != 0)
 824     {
 825       PsiVector::Iterator iter(*(this->array[c]));
 826       usint pos = iter.valueAfter(range.first).first;
 827       while(pos <= range.second)
 828       {
 829         bwt[pos - range.first] = c;
 830         pos = iter.selectNext();
 831       }
 832     }
 833   }
 834
 835   return bwt;
 836 }
 837
 838 usint
 839 RLCSA::countRuns() const
 840 {
 841   usint runs = 0;
 842   for(usint c = 0; c < CHARS; c++)
 843   {
 844     if(this->array[c] != 0)
 845     {
 846       PsiVector::Iterator iter(*(this->array[c]));
 847       runs += iter.countRuns();
 848     }
 849   }
 850
 851   return runs;
 852 }
 853
 854 //--------------------------------------------------------------------------
 855
 856 PsiVector::Iterator**
 857 RLCSA::getIterators() const
 858 {
 859   PsiVector::Iterator** iters = new PsiVector::Iterator*[CHARS];
 860   for(usint i = 0; i < CHARS; i++)
 861   {
 862     if(this->array[i] == 0) { iters[i] = 0; }
 863     else                    { iters[i] = new PsiVector::Iterator(*(this->array[i])); }
 864   }
 865   return iters;
 866 }
 867
 868 void
 869 RLCSA::deleteIterators(PsiVector::Iterator** iters) const
 870 {
 871   if(iters == 0) { return; }
 872   for(usint i = 0; i < CHARS; i++) { delete iters[i]; }
 873   delete[] iters;
 874 }
 875
 876 //--------------------------------------------------------------------------
 877
 878 usint
 879 RLCSA::reportSize(bool print) const
 880 {
 881   usint bytes = 0, temp = 0, bwt = 0;
 882
 883   for(usint c = 0; c < CHARS; c++)
 884   {
 885     if(this->array[c])
 886     {
 887       bytes += this->array[c]->reportSize();
 888       bwt += this->array[c]->getCompressedSize();
 889     }
 890   }
 891   bytes += sizeof(*this) + this->end_points->reportSize();
 892   if(print)
 893   {
 894     std::cout << "RLCSA:           " << (bytes / (double)MEGABYTE) << " MB" << std::endl;
 895     std::cout << "  BWT only:      " << (bwt / (double)MEGABYTE) << " MB" << std::endl;
 896   }
 897
 898   if(this->support_locate || this->support_display)
 899   {
 900     temp = this->sa_samples->reportSize();
 901     if(print) { std::cout << "SA samples:      " << (temp / (double)MEGABYTE) << " MB" << std::endl; }
 902     bytes += temp;
 903   }
 904
 905   if(print)
 906   {
 907     std::cout << "Total size:      " << (bytes / (double)MEGABYTE) << " MB" << std::endl;
 908     std::cout << std::endl;
 909   }
 910
 911   return bytes;
 912 }
 913
 914 void
 915 RLCSA::printInfo() const
 916 {
 917   double megabytes = this->data_size / (double)MEGABYTE;
 918
 919   std::cout << "Sequences:       " << this->number_of_sequences << std::endl;
 920   std::cout << "Original size:   " << megabytes << " MB" << std::endl;
 921   std::cout << "Block size:      " << (this->getBlockSize() * sizeof(usint)) << " bytes" << std::endl;
 922   if(this->support_locate || this->support_display)
 923   {
 924     std::cout << "Sample rate:     " << this->sample_rate << std::endl;
 925   }
 926   std::cout << std::endl;
 927 }
 928
 929 //--------------------------------------------------------------------------
 930
 931 RLEVector*
 932 RLCSA::buildPLCP(usint block_size) const
 933 {
 934   if(block_size < 2 * sizeof(usint) || block_size % sizeof(usint) != 0)
 935   {
 936     std::cerr << "PLCP: Block size must be a multiple of " << sizeof(usint) << " bytes!" << std::endl;
 937     return 0;
 938   }
 939
 940   PsiVector::Iterator** iters = this->getIterators();
 941
 942   RLEEncoder plcp(block_size);
 943   std::list<pair_type> matches;
 944   pair_type prev_range = EMPTY_PAIR;
 945   for(usint j = 0; j < this->number_of_sequences; j++)
 946   {
 947     // Encode the padding as a single run.
 948     pair_type seq_range = this->getSequenceRange(j);
 949     if(j > 0 && prev_range.second + 1 < seq_range.first)
 950     {
 951       this->encodePLCPRun(plcp, prev_range.second + 1, seq_range.first, seq_range.first - prev_range.second - 2);
 952     }
 953     prev_range = seq_range;
 954
 955     usint maximal = seq_range.first;
 956     usint x = this->sa_samples->inverseSA(seq_range.first), next_x;
 957
 958     // Invariant: x == inverseSA(i)
 959     for(usint i = seq_range.first; i <= seq_range.second; i++, x = next_x)
 960     {
 961       usint c = this->getCharacter(x);
 962
 963       // T[i,n] is lexicographically the first suffix beginning with c.
 964       if(x == this->index_ranges[c].first)
 965       {
 966         if(!matches.empty())
 967         {
 968           maximal = (*(matches.begin())).first;
 969           matches.clear();
 970         }
 971         this->encodePLCPRun(plcp, maximal, i + 1, i - maximal);
 972         maximal = i + 1;
 973         next_x = this->psiUnsafe(x, c, *(iters[c])) - this->number_of_sequences;
 974         continue;
 975       }
 976
 977       // Process the previous left matches we are still following.
 978       usint low = this->index_ranges[c].first + this->number_of_sequences;
 979       usint start, stop; start = stop = maximal;
 980       for(std::list<pair_type>::iterator iter = matches.begin(); iter != matches.end();)
 981       {
 982         pair_type match = *iter;
 983         if(match.second < low)  // These no longer match the current character.
 984         {
 985           if(match.first < start) { start = match.first; }  // Form a single run from then.
 986           iter = matches.erase(iter);
 987         }
 988         else
 989         {
 990           if(match.first < stop) { stop = match.first; }  // End of the run to be encoded.
 991           (*iter).second = this->psiUnsafe(match.second - this->number_of_sequences, c, *(iters[c]));
 992           ++iter;
 993         }
 994       }
 995       if(start < stop) { this->encodePLCPRun(plcp, start, stop, i - start); }
 996
 997       // If PLCP[i] is minimal, we add the left match to the list.
 998       // We add pairs of type (j, inverseSA(j) + n_of_s) as inverseSA is < 0 for end markers.
 999       usint next_y = this->psiUnsafe(x - 1, c, *(iters[c]));
1000       next_x = this->psiUnsafeNext(c, *(iters[c])) - this->number_of_sequences;
1001       if(next_y != next_x + this->number_of_sequences - 1)
1002       {
1003         matches.push_back(pair_type(maximal, next_y));
1004         maximal = i + 1;
1005       }
1006     }
1007
1008     if(!matches.empty())
1009     {
1010       maximal = (*(matches.begin())).first;
1011       matches.clear();
1012     }
1013     if(maximal <= seq_range.second)
1014     {
1015       this->encodePLCPRun(plcp, maximal, seq_range.second + 1, seq_range.second + 1 - maximal);
1016     }
1017   }
1018
1019   this->deleteIterators(iters);
1020
1021   return new RLEVector(plcp, 2 * (this->end_points->getSize() + 1));
1022 }
1023
1024 usint
1025 RLCSA::sampleLCP(usint sample_rate, pair_type*& sampled_values, bool report) const
1026 {
1027   if(sample_rate == 0)
1028   {
1029     sample_rate = this->data_size + 1;
1030   }
1031
1032   PsiVector::Iterator** iters = this->getIterators();
1033
1034   usint runs = this->countRuns();
1035   usint samples = 0, minimal_sum = 0, nonstrict_minimal_sum = 0, minimal_samples = 0;
1036   usint max_samples = runs + (this->data_size - runs) / sample_rate;
1037   sampled_values = new pair_type[max_samples];
1038
1039   std::list<Triple> matches;
1040   for(usint j = 0; j < this->number_of_sequences; j++)
1041   {
1042     pair_type seq_range = this->getSequenceRange(j);
1043     usint first_sample = samples; // First minimal sample of the current sequence.
1044     usint minimal = seq_range.second + 1; // Last minimal position.
1045     usint i, x = this->sa_samples->inverseSA(seq_range.first), next_x, prev_x = 0;
1046
1047     // Invariant: x == inverseSA(i)
1048     for(i = seq_range.first; i <= seq_range.second; i++, x = next_x)
1049     {
1050       usint c = this->getCharacter(x);
1051
1052       // T[i,n] is lexicographically the first suffix beginning with c.
1053       if(x == this->index_ranges[c].first)
1054       {
1055         if(!matches.empty())
1056         {
1057           for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end(); ++iter)
1058           {
1059             nonstrict_minimal_sum += i - (*iter).first;
1060           }
1061           matches.clear();
1062         }
1063
1064         // Minimal sample: SA[x] = i, LCP[x] = 0
1065         sampled_values[samples] = pair_type(x, 0); minimal = i;
1066         samples++; minimal_samples++;
1067         next_x = this->psiUnsafe(x, c, *(iters[c])) - this->number_of_sequences; prev_x = x;
1068         continue;
1069       }
1070
1071       // Process the previous left matches we are still following.
1072       usint low = this->index_ranges[c].first + this->number_of_sequences;
1073       pair_type potential_sample(this->data_size, this->data_size);
1074       for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end();)
1075       {
1076         Triple match = *iter;
1077         if(match.second < low)  // These no longer match the current character.
1078         {
1079           if(potential_sample.first != this->data_size) { nonstrict_minimal_sum += potential_sample.second; }
1080
1081           // Potential minimal sample: SA[match.third] = match.first, LCP[match.third] = i - match.first
1082           potential_sample = pair_type(match.third, i - match.first);
1083           iter = matches.erase(iter);
1084         }
1085         else
1086         {
1087           (*iter).second = this->psiUnsafe(match.second - this->number_of_sequences, c, *(iters[c]));
1088           ++iter;
1089         }
1090       }
1091       if(potential_sample.first != this->data_size)
1092       {
1093         // Last potential sample is minimal.
1094         sampled_values[samples] = potential_sample;
1095         samples++; minimal_sum += potential_sample.second; minimal_samples++;
1096       }
1097
1098       // If PLCP[i] is minimal, we add the left match to the list.
1099       // We add triples of type (j, inverseSA(j) + n_of_s, x) as inverseSA is < 0 for end markers.
1100       // x is the original minimal position.
1101       usint next_y = this->psiUnsafe(x - 1, c, *(iters[c]));
1102       next_x = this->psiUnsafeNext(c, *(iters[c]));
1103       if(next_y != next_x - 1 || next_x < this->number_of_sequences)
1104       {
1105         matches.push_back(Triple(i, next_y, x)); minimal = i;
1106       }
1107       next_x -= this->number_of_sequences; prev_x = x;
1108     }
1109
1110     // Are we still following something?
1111     if(!matches.empty())
1112     {
1113       pair_type potential_sample(this->data_size, this->data_size);
1114       for(std::list<Triple>::iterator iter = matches.begin(); iter != matches.end(); ++iter)
1115       {
1116         Triple match = *iter;
1117         if(potential_sample.first != this->data_size) { nonstrict_minimal_sum += potential_sample.second; }
1118
1119         // Potential minimal sample: SA[match.third] = match.first, LCP[match.third] = i - match.first
1120         potential_sample = pair_type(match.third, i - match.first);
1121       }
1122       matches.clear();
1123       if(potential_sample.first != this->data_size)
1124       {
1125         // Last potential sample is minimal.
1126         sampled_values[samples] = potential_sample;
1127         samples++; minimal_sum += potential_sample.second; minimal_samples++;
1128       }
1129     }
1130
1131     // Add the non-minimal samples.
1132     if(sample_rate <= this->data_size)
1133     {
1134       usint last_sample = samples - 1;
1135       i = seq_range.first; x = this->sa_samples->inverseSA(seq_range.first);
1136       for(usint current_sample = first_sample; current_sample <= last_sample; current_sample++)
1137       {
1138         // Find the next minimal sample and add nonminimal samples if needed.
1139         pair_type first_nonminimal(i + sample_rate - 1, samples);
1140         pair_type next_nonminimal = first_nonminimal;
1141         while(x != sampled_values[current_sample].first)
1142         {
1143           if(i == next_nonminimal.first)
1144           {
1145             sampled_values[samples] = pair_type(x, 0); samples++;
1146             next_nonminimal.first += sample_rate; next_nonminimal.second++;
1147           }
1148           i++; x = this->psi(x, iters) - this->number_of_sequences;
1149         }
1150
1151         // Reduce the nonminimal samples to the current minimal sample.
1152         for(next_nonminimal = first_nonminimal; next_nonminimal.second < samples; next_nonminimal.first += sample_rate, next_nonminimal.second++)
1153         {
1154           sampled_values[next_nonminimal.second].second = sampled_values[current_sample].second + i - next_nonminimal.first;
1155         }
1156         i++; x = this->psi(x, iters) - this->number_of_sequences;
1157       }
1158     }
1159   }
1160
1161   std::sort(sampled_values, sampled_values + samples);
1162
1163   if(report)
1164   {
1165     std::cout << "Samples: " << samples << " (total) / " << minimal_samples << " (minimal)" << std::endl;
1166     std::cout << "Upper bounds: " << max_samples << " (total) / " << runs << " (minimal)" << std::endl;
1167     std::cout << "Sum of minimal samples: " << (minimal_sum + nonstrict_minimal_sum) << " (total) / " << minimal_sum << " (strict)" << std::endl;
1168     std::cout << std::endl;
1169   }
1170
1171   this->deleteIterators(iters);
1172
1173   return samples;
1174 }
1175
1176 usint
1177 RLCSA::lcp(usint index, const LCPSamples& lcp_samples) const
1178 {
1179   if(index >= this->data_size) { return this->data_size; }
1180
1181   usint offset = 0;
1182   index += this->number_of_sequences;
1183   while(true)
1184   {
1185     pair_type next_sample = lcp_samples.getFirstSampleAfter(index - this->number_of_sequences);
1186     next_sample.first += this->number_of_sequences;
1187     if(next_sample.first == index)
1188     {
1189       return lcp_samples.getSample(next_sample.second) + offset;
1190     }
1191     index = this->psi(index - this->number_of_sequences);
1192     offset++;
1193   }
1194 }
1195
1196 usint
1197 RLCSA::lcp(usint index, const LCPSamples& lcp_samples, const RLEVector& plcp) const
1198 {
1199   if(index >= this->data_size) { return this->data_size; }
1200
1201   usint offset = 0;
1202   while(true)
1203   {
1204     usint c = this->getCharacter(index);
1205     PsiVector::Iterator iter(*(this->array[c]));
1206
1207     if(index == this->index_ranges[c].first)
1208     {
1209       return lcp_samples.getSampleAtPosition(index) + offset;
1210     }
1211     usint next_y = this->psiUnsafe(index - 1, c, iter);
1212     usint next_x = this->psiUnsafeNext(c, iter);
1213     if(next_y != next_x - 1 || next_x < this->number_of_sequences)
1214     {
1215       pair_type sample = lcp_samples.getFirstSampleAfter(index);
1216       if(sample.first == index)
1217       {
1218         return lcp_samples.getSample(sample.second) + offset;
1219       }
1220     }
1221
1222     pair_type next_sample = this->sa_samples->getFirstSampleAfter(index);
1223     if(next_sample.first == index)
1224     {
1225       index = this->sa_samples->getSample(next_sample.second) - offset;
1226       RLEVector::Iterator plcp_iter(plcp);
1227       return plcp_iter.select(index) - 2 * index;
1228     }
1229     index = next_x - this->number_of_sequences;
1230     offset++;
1231   }
1232 }
1233
1234 usint
1235 RLCSA::lcpDirect(usint index) const
1236 {
1237   if(index >= this->data_size || index == 0)
1238   {
1239     return 0;
1240   }
1241
1242   usint match = index - 1;
1243   usint value = 0;
1244
1245   usint c = this->getCharacter(index);
1246   usint low = this->index_ranges[c].first;
1247   while(match >= low)
1248   {
1249     PsiVector::Iterator iter(*(this->array[c]));
1250     match = this->psiUnsafe(match, c, iter);
1251     index = this->psiUnsafe(index, c, iter);
1252     value++;
1253
1254     if(match < this->number_of_sequences || index < this->number_of_sequences)
1255     {
1256       break;
1257     }
1258     match -= this->number_of_sequences;
1259     index -= this->number_of_sequences;
1260
1261     c = this->getCharacter(index);
1262     low = this->index_ranges[c].first;
1263   }
1264
1265   return value;
1266 }
1267
1268 //--------------------------------------------------------------------------
1269
1270 void
1271 RLCSA::mergeEndPoints(RLCSA& index, RLCSA& increment)
1272 {
1273   DeltaEncoder* endings = new DeltaEncoder(RLCSA::ENDPOINT_BLOCK_SIZE);
1274
1275   DeltaVector::Iterator index_iter(*(index.end_points));
1276   DeltaVector::Iterator increment_iter(*(increment.end_points));
1277
1278   endings->setBit(index_iter.select(0));
1279   for(usint i = 1; i < index.number_of_sequences; i++)
1280   {
1281     endings->setBit(index_iter.selectNext());
1282   }
1283   usint sum = index.end_points->getSize();
1284   delete index.end_points; index.end_points = 0;
1285
1286   endings->setBit(sum + increment_iter.select(0));
1287   for(usint i = 1; i < increment.number_of_sequences; i++)
1288   {
1289     endings->setBit(sum + increment_iter.selectNext());
1290   }
1291   sum += increment.end_points->getSize();
1292   delete increment.end_points; increment.end_points = 0;
1293
1294   this->end_points = new DeltaVector(*endings, sum);
1295   delete endings;
1296 }
1297
1298
1299 void
1300 RLCSA::mergeSamples(RLCSA& index, RLCSA& increment, usint* positions)
1301 {
1302   if(this->support_locate || this->support_display)
1303   {
1304     positions += increment.number_of_sequences;
1305     this->sa_samples = new SASamples(*(index.sa_samples), *(increment.sa_samples), positions, increment.data_size, this->number_of_sequences);
1306   }
1307 }
1308
1309 //--------------------------------------------------------------------------
1310
1311 void
1312 RLCSA::buildCharIndexes(usint* distribution)
1313 {
1314   this->data_size = buildRanges(distribution, this->index_ranges);
1315
1316   usint i = 0, c = 0;
1317   for(; c < CHARS; c++)
1318   {
1319     if(!isEmpty(this->index_ranges[c]))
1320     {
1321       this->text_chars[i] = c;
1322       i++;
1323     }
1324   }
1325   this->chars = i;
1326
1327   this->index_rate = std::max((this->data_size + CHARS - 1) / CHARS, (usint)1);
1328   usint current = 0;
1329
1330   for(c = 0, i = 0; c < this->chars; c++)
1331   {
1332     pair_type range = this->index_ranges[this->text_chars[c]];
1333     while(current <= range.second)
1334     {
1335       this->index_pointers[i] = c;
1336       current += this->index_rate;
1337       i++;
1338     }
1339   }
1340 }
1341
1342 usint
1343 buildRanges(usint* distribution, pair_type* index_ranges)
1344 {
1345   if(distribution == 0 || index_ranges == 0) { return 0; }
1346
1347   usint sum = 0;
1348   for(usint c = 0; c < CHARS; c++)
1349   {
1350     if(distribution[c] == 0)
1351     {
1352       if(sum == 0) { index_ranges[c] = EMPTY_PAIR; }
1353       else         { index_ranges[c] = pair_type(sum, sum - 1); }
1354     }
1355     else
1356     {
1357       index_ranges[c].first = sum;
1358       sum += distribution[c];
1359       index_ranges[c].second = sum - 1;
1360     }
1361   }
1362
1363   return sum;
1364 }
1365
1366
1367 } // namespace CSA