incbwt/rlcsa.cpp

   1 #include <algorithm>
   2 #include <cstdlib>
   3 #include <ctime>
   4 #include <fstream>
   5 #include <iostream>
   6
   7 #include "rlcsa.h"
   8 #include "misc/utils.h"
   9 #include "qsufsort/qsufsort.h"
  10
  11
  12 namespace CSA
  13 {
  14
  15
  16 RLCSA::RLCSA(const std::string& base_name, bool print) :
  17   ok(false),
  18   sa_samples(0),
  19   end_points(0)
  20 {
  21   for(usint c = 0; c < CHARS; c++) { this->array[c] = 0; }
  22
  23   std::string array_name = base_name + ARRAY_EXTENSION;
  24   std::ifstream array_file(array_name.c_str(), std::ios_base::binary);
  25   if(!array_file)
  26   {
  27     std::cerr << "RLCSA: Error opening Psi array file!" << std::endl;
  28     return;
  29   }
  30
  31   usint distribution[CHARS];
  32   array_file.read((char*)distribution, CHARS * sizeof(usint));
  33   this->buildCharIndexes(distribution);
  34
  35   Parameters parameters;
  36   parameters.read(base_name + PARAMETERS_EXTENSION);
  37   for(usint c = 0; c < CHARS; c++)
  38   {
  39     if(!isEmpty(this->index_ranges[c])) { this->array[c] = new RLEVector(array_file); }
  40   }
  41
  42   this->end_points = new DeltaVector(array_file);
  43   this->number_of_sequences = this->end_points->getNumberOfItems();
  44
  45   array_file.read((char*)&(this->sample_rate), sizeof(this->sample_rate));
  46   array_file.close();
  47
  48   this->support_locate = parameters.get(SUPPORT_LOCATE);
  49   this->support_display = parameters.get(SUPPORT_DISPLAY);
  50
  51   if(this->support_locate || this->support_display)
  52   {
  53     std::string sa_sample_name = base_name + SA_SAMPLES_EXTENSION;
  54     std::ifstream sa_sample_file(sa_sample_name.c_str(), std::ios_base::binary);
  55     if(!sa_sample_file)
  56     {
  57       std::cerr << "RLCSA: Error opening suffix array sample file!" << std::endl;
  58       return;
  59     }
  60     this->sa_samples = new SASamples(sa_sample_file, this->sample_rate);
  61     sa_sample_file.close();
  62   }
  63
  64   if(print) { parameters.print(); }
  65
  66   this->ok = true;
  67 }
  68
  69 RLCSA::RLCSA(uchar* data, usint bytes, usint block_size, usint sa_sample_rate, bool multiple_sequences, bool delete_data) :
  70   ok(false),
  71   sa_samples(0),
  72   sample_rate(sa_sample_rate),
  73   end_points(0)
  74 {
  75   for(usint c = 0; c < CHARS; c++) { this->array[c] = 0; }
  76
  77   if(!data || bytes == 0)
  78   {
  79     std::cerr << "RLCSA: No input data given!" << std::endl;
  80     if(delete_data) { delete[] data; }
  81     return;
  82   }
  83   if(block_size < 2 * sizeof(usint) || block_size % sizeof(usint) != 0)
  84   {
  85     std::cerr << "RLCSA: Block size must be a multiple of " << sizeof(usint) << " bytes!" << std::endl;
  86     if(delete_data) { delete[] data; }
  87     return;
  88   }
  89
  90
  91   // Do we store SA samples?
  92   if(this->sample_rate > 0)
  93   {
  94     this->support_locate = this->support_display = true;
  95   }
  96   else
  97   {
  98     this->support_locate = this->support_display = false;
  99     this->sample_rate = 1;
 100   }
 101
 102
 103   // Determine the number of sequences and mark their end points.
 104   if(multiple_sequences)
 105   {
 106     DeltaEncoder endings(RLCSA::ENDPOINT_BLOCK_SIZE);
 107     this->number_of_sequences = 0;
 108     usint marker = 0;
 109     usint padding = 0, chars_encountered = 0;
 110
 111     for(usint i = 0; i < bytes; i++)
 112     {
 113       if(data[i] == 0)
 114       {
 115         if(i == marker) { break; }  // Empty sequence.
 116         this->number_of_sequences++;
 117         marker = i + 1;
 118
 119         usint pos = chars_encountered + padding - 1;
 120         endings.setBit(pos);
 121         padding = ((pos + this->sample_rate) / this->sample_rate) * this->sample_rate - chars_encountered;
 122       }
 123       else
 124       {
 125         chars_encountered++;
 126       }
 127     }
 128
 129     if(this->number_of_sequences == 0 || marker != bytes)
 130     {
 131       std::cerr << "RLCSA: Collection must consist of 0-terminated nonempty sequences!" << std::endl;
 132       if(delete_data) { delete[] data; }
 133       return;
 134     }
 135     this->end_points = new DeltaVector(endings, chars_encountered + padding);
 136   }
 137   else
 138   {
 139     this->number_of_sequences = 1;
 140     DeltaEncoder endings(RLCSA::ENDPOINT_BLOCK_SIZE, RLCSA::ENDPOINT_BLOCK_SIZE);
 141     usint pos = bytes - 1;
 142     endings.setBit(pos);
 143     usint padding = ((pos + this->sample_rate) / this->sample_rate) * this->sample_rate - bytes;
 144     this->end_points = new DeltaVector(endings, bytes + padding);
 145   }
 146
 147
 148   // Build character tables etc.
 149   usint distribution[CHARS];
 150   sint low = CHARS, high = 0;
 151   for(usint c = 0; c < CHARS; c++) { distribution[c] = 0; }
 152   for(usint i = 0; i < bytes; i++)
 153   {
 154     if(data[i] < low) { low = data[i]; }
 155     if(data[i] > high) { high = data[i]; }
 156     distribution[(usint)data[i]]++;
 157   }
 158   if(multiple_sequences)
 159   {
 160     distribution[0] = 0;
 161     low = 0;
 162     high = high + this->number_of_sequences;
 163   }
 164   this->buildCharIndexes(distribution);
 165
 166
 167   // Build suffix array.
 168   sint* inverse = new sint[bytes + 1];
 169   if(multiple_sequences)
 170   {
 171     sint zeros = 0;
 172     for(usint i = 0; i < bytes; i++)
 173     {
 174       if(data[i] == 0)
 175       {
 176         inverse[i] = zeros;
 177         zeros++;
 178       }
 179       else
 180       {
 181         inverse[i] = (sint)data[i] + this->number_of_sequences;
 182       }
 183     }
 184   }
 185   else
 186   {
 187     for(usint i = 0; i < bytes; i++) { inverse[i] = (sint)data[i]; }
 188   }
 189   if(delete_data) { delete[] data; }
 190   sint* sa = new sint[bytes + 1];
 191   suffixsort(inverse, sa, bytes, high + 1, low);
 192
 193
 194   // Sample SA.
 195   usint incr = (multiple_sequences ? this->number_of_sequences + 1 : 1);  // No e_of_s markers in SA.
 196   if(this->support_locate || this->support_display)
 197   {
 198     if(multiple_sequences)
 199     {
 200       std::cout << "We shouldn't be here!" << std::endl;
 201       // Real SA starts at sa + incr.
 202       // for each sequence
 203       //   find starting position
 204       //   sample relative to starting position
 205       // sort samples to SA order
 206     }
 207     else
 208     {
 209       this->sa_samples = new SASamples((usint*)(sa + incr), this->data_size, this->sample_rate);
 210     }
 211   }
 212
 213
 214   // Build Psi.
 215   for(usint i = 0; i <= bytes; i++)
 216   {
 217     sa[i] = inverse[(sa[i] + 1) % (bytes + 1)];
 218   }
 219   delete[] inverse;
 220
 221
 222   // Build RLCSA.
 223   usint decr = (multiple_sequences ? 1 : 0);  // Ignore the last e_of_s marker if multiple sequences.
 224   for(usint c = 0; c < CHARS; c++)
 225   {
 226     if(distribution[c] == 0) { this->array[c] = 0; continue; }
 227
 228     usint* curr = (usint*)(sa + index_ranges[c].first + incr);
 229     usint* limit = (usint*)(sa + index_ranges[c].second + incr + 1);
 230     RLEEncoder encoder(block_size);
 231     pair_type run(*curr, 1);
 232     curr++;
 233
 234     for(; curr < limit; curr++)
 235     {
 236       if(*curr == run.first + run.second) { run.second++; }
 237       else
 238       {
 239         encoder.setRun(run.first - decr, run.second);
 240         run = pair_type(*curr, 1);
 241       }
 242     }
 243     encoder.setRun(run.first - decr, run.second);
 244
 245     this->array[c] = new RLEVector(encoder, this->data_size + incr - decr);
 246   }
 247   delete[] sa;
 248
 249
 250   this->ok = true;
 251 }
 252
 253 RLCSA::RLCSA(RLCSA& index, RLCSA& increment, usint* positions, usint block_size) :
 254   ok(false),
 255   sa_samples(0),
 256   end_points(0)
 257 {
 258   for(usint c = 0; c < CHARS; c++) { this->array[c] = 0; }
 259
 260   if(!index.isOk() || !increment.isOk())
 261   {
 262     return; // Fail silently. Actual error has already been reported.
 263   }
 264   if(positions == 0)
 265   {
 266     std::cerr << "RLCSA: Positions for insertions not available!" << std::endl;
 267     return;
 268   }
 269   if(index.sample_rate != increment.sample_rate)
 270   {
 271     std::cerr << "RLCSA: Cannot combine indexes with different sample rates!" << std::endl;
 272     std::cout << "Index: " << index.sample_rate << ", increment: " << increment.sample_rate << std::endl;
 273     return;
 274   }
 275
 276   if(index.sa_samples == 0 || increment.sa_samples == 0)
 277   {
 278     this->support_locate = this->support_display = false;
 279   }
 280   else
 281   {
 282     this->support_locate = this->support_display = true;
 283   }
 284
 285
 286   // Build character tables etc.
 287   usint distribution[CHARS];
 288   for(usint c = 0; c < CHARS; c++)
 289   {
 290     distribution[c] = length(index.index_ranges[c]) + length(increment.index_ranges[c]);
 291   }
 292   this->buildCharIndexes(distribution);
 293   this->sample_rate = index.sample_rate;
 294
 295
 296   // Combine end points of sequences.
 297   this->number_of_sequences = index.number_of_sequences + increment.number_of_sequences;
 298   DeltaEncoder* endings = new DeltaEncoder(RLCSA::ENDPOINT_BLOCK_SIZE);
 299
 300   endings->setBit(index.end_points->select(0));
 301   for(usint i = 1; i < index.number_of_sequences; i++)
 302   {
 303     endings->setBit(index.end_points->selectNext());
 304   }
 305   usint sum = index.end_points->getSize();
 306   delete index.end_points; index.end_points = 0;
 307
 308   endings->setBit(sum + increment.end_points->select(0));
 309   for(usint i = 1; i < increment.number_of_sequences; i++)
 310   {
 311     endings->setBit(sum + increment.end_points->selectNext());
 312   }
 313   sum += increment.end_points->getSize();
 314   delete increment.end_points; increment.end_points = 0;
 315
 316   this->end_points = new DeltaVector(*endings, sum);
 317   delete endings;
 318
 319
 320   // Combine Psi.
 321   usint psi_size = this->data_size + this->number_of_sequences;
 322   for(usint c = 0; c < CHARS; c++)
 323   {
 324     if(distribution[c] == 0) { this->array[c] = 0; continue; }
 325     this->array[c] = mergeVectors(index.array[c], increment.array[c], positions, increment.data_size + increment.number_of_sequences, psi_size, block_size);
 326     index.array[c] = 0;
 327     increment.array[c] = 0;
 328
 329     if(this->array[c] == 0)
 330     {
 331       std::cerr << "RLCSA: Merge failed for vectors " << c << "!" << std::endl;
 332       return;
 333     }
 334   }
 335
 336
 337   // Combine suffix array samples.
 338   if(this->support_locate || this->support_display)
 339   {
 340     positions += increment.number_of_sequences;
 341     for(usint i = 0; i < increment.data_size; i++)
 342     {
 343       positions[i] -= this->number_of_sequences;
 344     }
 345     this->sa_samples = new SASamples(*(index.sa_samples), *(increment.sa_samples), positions, increment.data_size);
 346   }
 347
 348   this->ok = true;
 349 }
 350
 351 RLCSA::~RLCSA()
 352 {
 353   for(usint c = 0; c < CHARS; c++) { delete this->array[c]; }
 354   delete this->sa_samples;
 355   delete this->end_points;
 356 }
 357
 358 void
 359 RLCSA::writeTo(const std::string& base_name)
 360 {
 361   std::string array_name = base_name + ARRAY_EXTENSION;
 362   std::ofstream array_file(array_name.c_str(), std::ios_base::binary);
 363   if(!array_file)
 364   {
 365     std::cerr << "RLCSA: Error creating Psi array file!" << std::endl;
 366     return;
 367   }
 368
 369   usint distribution[CHARS];
 370   for(usint c = 0; c < CHARS; c++)
 371   {
 372     distribution[c] = length(this->index_ranges[c]);
 373   }
 374   array_file.write((char*)distribution, CHARS * sizeof(usint));
 375   for(usint c = 0; c < CHARS; c++)
 376   {
 377     if(this->array[c] != 0)
 378     {
 379       this->array[c]->writeTo(array_file);
 380     }
 381   }
 382
 383   this->end_points->writeTo(array_file);
 384   array_file.write((char*)&(this->sample_rate), sizeof(this->sample_rate));
 385   array_file.close();
 386
 387   if(this->support_locate || this->support_display)
 388   {
 389     std::string sa_sample_name = base_name + SA_SAMPLES_EXTENSION;
 390     std::ofstream sa_sample_file(sa_sample_name.c_str(), std::ios_base::binary);
 391     if(!sa_sample_file)
 392     {
 393       std::cerr << "RLCSA: Error creating suffix array sample file!" << std::endl;
 394       return;
 395     }
 396
 397     this->sa_samples->writeTo(sa_sample_file);
 398     sa_sample_file.close();
 399   }
 400 }
 401
 402 bool
 403 RLCSA::isOk()
 404 {
 405   return this->ok;
 406 }
 407
 408 //--------------------------------------------------------------------------
 409
 410 pair_type
 411 RLCSA::count(const std::string& pattern)
 412 {
 413   if(pattern.length() == 0) { return pair_type(0, this->data_size - 1); }
 414
 415   pair_type index_range = this->index_ranges[(usint)*(pattern.rbegin())];
 416   if(isEmpty(index_range)) { return index_range; }
 417   index_range.first += this->number_of_sequences;
 418   index_range.second += this->number_of_sequences;
 419
 420   for(std::string::const_reverse_iterator iter = ++pattern.rbegin(); iter != pattern.rend(); iter++)
 421   {
 422     RLEVector* vector = this->array[(usint)*iter];
 423     usint start = this->index_ranges[(usint)*iter].first;
 424
 425     index_range.first = start + vector->rank(index_range.first, true) - 1 + this->number_of_sequences;
 426     index_range.second = start + vector->rank(index_range.second) - 1 + this->number_of_sequences;
 427
 428     if(isEmpty(index_range)) { return index_range; }
 429   }
 430
 431   // Suffix array indexes are 0-based.
 432   index_range.first -= this->number_of_sequences;
 433   index_range.second -= this->number_of_sequences;
 434
 435   return index_range;
 436 }
 437
 438 void
 439 RLCSA::reportPositions(uchar* data, usint length, usint* positions)
 440 {
 441   if(data == 0 || length == 0 || positions == 0) { return; }
 442
 443   usint current = this->number_of_sequences - 1;
 444   positions[length] = current; // "immediately after current"
 445   for(sint i = (sint)(length - 1); i >= 0; i--)
 446   {
 447 //     positions[i] = current; // "immediately after current"
 448     usint c = (usint)data[i];
 449     if(array[c])
 450     {
 451       current = this->index_ranges[c].first + this->array[c]->rank(current) - 1 + this->number_of_sequences;
 452     }
 453     else
 454     {
 455       if(c < this->text_chars[0]) // No previous characters either.
 456       {
 457         current = this->number_of_sequences - 1;
 458       }
 459       else
 460       {
 461         current = this->index_ranges[c].first - 1 + this->number_of_sequences;
 462       }
 463     }
 464     positions[i] = current; // "immediately after current"
 465   }
 466 }
 467
 468 //--------------------------------------------------------------------------
 469
 470 LocateItem*
 471 RLCSA::locate(pair_type range)
 472 {
 473   if(!(this->support_locate) || isEmpty(range) || range.second >= this->data_size) { return 0; }
 474
 475   LocateItem* data = new LocateItem[range.second + 1 - range.first];
 476   if(!data) { return 0; }
 477   this->locateUnsafe(range, data);
 478
 479   return data;
 480 }
 481
 482 LocateItem*
 483 RLCSA::locate(pair_type range, LocateItem* data)
 484 {
 485   if(!(this->support_locate) || isEmpty(range) || range.second >= this->data_size || data == 0) { return 0; }
 486   this->locateUnsafe(range, data);
 487   return data;
 488 }
 489
 490 void
 491 RLCSA::locateUnsafe(pair_type range, LocateItem* data)
 492 {
 493   usint items = range.second + 1 - range.first;
 494   for(usint i = 0, j = range.first; i < items; i++, j++)
 495   {
 496     data[i].value = j + this->number_of_sequences;
 497     data[i].offset = 0;
 498     data[i].found = false;
 499   }
 500
 501   bool found = false;
 502   while(!found)
 503   {
 504     found = true;
 505     pair_type run = EMPTY_PAIR;
 506     for(usint i = 0; i < items; i++)
 507     {
 508       if(data[i].found)
 509       {
 510         continue; // The run might continue after this.
 511       }
 512       else if(isEmpty(run))
 513       {
 514         run = pair_type(i, i);
 515       }
 516       else if(data[i].value - data[run.first].value == i - run.first)
 517       {
 518         run.second = i;
 519       }
 520       else
 521       {
 522         found &= this->processRun(run, data);
 523         run = pair_type(i, i);
 524       }
 525     }
 526     if(!isEmpty(run)) { found &= this->processRun(run, data); }
 527   }
 528 }
 529
 530 bool
 531 RLCSA::processRun(pair_type run, LocateItem* data)
 532 {
 533   bool found = true;
 534   usint run_start = 0, run_left = 0;
 535   pair_type next_sample = pair_type(0, 0);
 536
 537   for(usint i = run.first; i <= run.second; i++)
 538   {
 539     if(data[i].found)
 540     {
 541       if(run_left > 0) { run_left--; }
 542       continue;
 543     }
 544     if(data[i].value < this->number_of_sequences) // Implicit sample here.
 545     {
 546       data[i].value = this->end_points->select(data[i].value) + 1 - data[i].offset;
 547       data[i].found = true;
 548       if(run_left > 0) { run_left--; }
 549       continue;
 550     }
 551     if(next_sample.first < data[i].value) // Need another sample.
 552     {
 553       next_sample = this->sa_samples->getFirstSampleAfter(data[i].value - this->number_of_sequences);
 554       next_sample.first += this->number_of_sequences;
 555     }
 556     if(data[i].value < next_sample.first) // No sample found for current position.
 557     {
 558       if(run_left > 0)
 559       {
 560         data[i].value = data[run_start].value + i - run_start;
 561         run_left--;
 562       }
 563       else
 564       {
 565         pair_type value = this->psi(data[i].value - this->number_of_sequences, run.second - i);
 566         data[i].value = value.first;
 567         run_left = value.second;
 568         run_start = i;
 569       }
 570       data[i].offset++;
 571       found = false;
 572     }
 573     else  // Sampled position found.
 574     {
 575       data[i].value = this->sa_samples->getSample(next_sample.second) - data[i].offset;
 576       data[i].found = true;
 577       if(run_left > 0) { run_left--; }
 578     }
 579   }
 580   return found;
 581 }
 582
 583 usint
 584 RLCSA::locate(usint index)
 585 {
 586   if(!(this->support_locate) || index >= this->data_size) { return this->data_size; }
 587
 588   usint offset = 0;
 589   index += this->number_of_sequences;
 590   while(true)
 591   {
 592     if(index < this->number_of_sequences) // Implicit sample here
 593     {
 594       return this->end_points->select(index) + 1 - offset;
 595     }
 596     pair_type next_sample = this->sa_samples->getFirstSampleAfter(index - this->number_of_sequences);
 597     next_sample.first += this->number_of_sequences;
 598     if(next_sample.first == index)
 599     {
 600       return this->sa_samples->getSample(next_sample.second) - offset;
 601     }
 602     index = this->psi(index - this->number_of_sequences);
 603     offset++;
 604   }
 605 }
 606
 607 //--------------------------------------------------------------------------
 608
 609 uchar*
 610 RLCSA::display(usint sequence, pair_type range)
 611 {
 612   if(!(this->support_display) || isEmpty(range)) { return 0; }
 613
 614   pair_type seq_range = this->getSequence(sequence);
 615   if(isEmpty(seq_range)) { return 0; }
 616
 617   range.first += seq_range.first; range.second += seq_range.first;
 618   if(range.second > seq_range.second) { return 0; }
 619
 620   uchar* data = new uchar[range.second + 1 - range.first];
 621   if(!data) { return 0; }
 622   this->displayUnsafe(range, data);
 623
 624   return data;
 625 }
 626
 627 uchar*
 628 RLCSA::display(usint sequence, pair_type range, uchar* data)
 629 {
 630   if(!(this->support_display) || isEmpty(range) || data == 0) { return 0; }
 631
 632   pair_type seq_range = this->getSequence(sequence);
 633   if(isEmpty(seq_range)) { return 0; }
 634
 635   range.first += seq_range.first; range.second += seq_range.first;
 636   if(range.second > seq_range.second) { return 0; }
 637
 638   this->displayUnsafe(range, data);
 639   return data;
 640 }
 641
 642 void
 643 RLCSA::displayUnsafe(pair_type range, uchar* data)
 644 {
 645   usint i = range.first - range.first % this->sa_samples->getSampleRate();
 646
 647   usint pos = this->sa_samples->inverseSA(i);
 648
 649   for(; i < range.first; i++)
 650   {
 651     pos = this->psi(pos) - this->number_of_sequences;
 652   }
 653   for(; i <= range.second; i++)
 654   {
 655     data[i - range.first] = this->getCharacter(pos);
 656     pos = this->psi(pos) - this->number_of_sequences;
 657   }
 658 }
 659
 660 //--------------------------------------------------------------------------
 661
 662 void
 663 RLCSA::decompressInto(std::ofstream& psi_file)
 664 {
 665   for(usint c = 0; c < CHARS; c++)
 666   {
 667     if(!(this->array[c])) { continue; }
 668     usint value = this->array[c]->select(0);
 669     psi_file.write((char*)&(value), sizeof(value));
 670     while(this->array[c]->hasNext())
 671     {
 672       value = this->array[c]->selectNext();
 673       psi_file.write((char*)&value, sizeof(value));
 674     }
 675   }
 676 }
 677
 678 uchar*
 679 RLCSA::readBWT()
 680 {
 681   usint n = this->data_size + this->number_of_sequences;
 682
 683   uchar* bwt = new uchar[n];
 684   memset(bwt, 0, n);
 685
 686   for(usint c = 0; c < CHARS; c++)
 687   {
 688     RLEVector* vector = this->array[c];
 689     if(vector != 0)
 690     {
 691       bwt[vector->select(0)] = c;
 692       while(vector->hasNext()) { bwt[vector->selectNext()] = c; }
 693     }
 694   }
 695
 696   return bwt;
 697 }
 698
 699 //--------------------------------------------------------------------------
 700
 701 usint
 702 RLCSA::psi(usint index)
 703 {
 704   if(index >= this->data_size)
 705   {
 706     return this->data_size + this->number_of_sequences;
 707   }
 708
 709   usint c = this->getCharacter(index);
 710   return this->array[c]->select(index - this->index_ranges[c].first);
 711 }
 712
 713 pair_type
 714 RLCSA::psi(usint index, usint max_length)
 715 {
 716   if(index >= this->data_size)
 717   {
 718     return pair_type(this->data_size + this->number_of_sequences, 0);
 719   }
 720
 721   usint c = this->getCharacter(index);
 722   return this->array[c]->selectRun(index - this->index_ranges[c].first, max_length);
 723 }
 724
 725 //--------------------------------------------------------------------------
 726
 727 pair_type
 728 RLCSA::getSequence(usint number)
 729 {
 730   if(number >= this->number_of_sequences) { return EMPTY_PAIR; }
 731
 732   pair_type result;
 733   if(number == 0)
 734   {
 735     result.first = 0;
 736     result.second = this->end_points->select(number);
 737   }
 738   else
 739   {
 740     if(this->sa_samples == 0) { return EMPTY_PAIR; }
 741     usint d = this->sa_samples->getSampleRate();
 742     result.first = d * ((this->end_points->select(number - 1) / d) + 1);
 743     result.second = this->end_points->selectNext();
 744   }
 745
 746   return result;
 747 }
 748
 749 usint
 750 RLCSA::reportSize(bool print)
 751 {
 752   usint bytes = 0, temp = 0;
 753
 754   for(usint c = 0; c < CHARS; c++)
 755   {
 756     if(this->array[c]) { temp += this->array[c]->reportSize(); }
 757   }
 758   if(print) { std::cout << "Psi array:       " << temp << std::endl; }
 759   bytes += temp;
 760
 761   if(this->support_locate || this->support_display)
 762   {
 763     temp = this->sa_samples->reportSize();
 764     if(print) { std::cout << "SA samples:      " << temp << std::endl; }
 765     bytes += temp;
 766   }
 767
 768   temp = sizeof(*this) + this->end_points->reportSize();
 769   if(print) { std::cout << "RLCSA overhead:  " << temp << std::endl; }
 770   bytes += temp;
 771
 772   if(print)
 773   {
 774     std::cout << "Total size:      " << bytes << std::endl;
 775     std::cout << std::endl;
 776   }
 777
 778   return bytes;
 779 }
 780
 781 //--------------------------------------------------------------------------
 782
 783 void
 784 RLCSA::buildCharIndexes(usint* distribution)
 785 {
 786   this->data_size = buildRanges(distribution, this->index_ranges);
 787
 788   usint i = 0, c = 0;
 789   for(; c < CHARS; c++)
 790   {
 791     if(!isEmpty(this->index_ranges[c]))
 792     {
 793       this->text_chars[i] = c;
 794       i++;
 795     }
 796   }
 797   this->chars = i;
 798
 799   this->index_rate = std::max((this->data_size + CHARS - 1) / CHARS, (usint)1);
 800   usint current = 0;
 801
 802   for(c = 0, i = 0; c < this->chars; c++)
 803   {
 804     pair_type range = this->index_ranges[this->text_chars[c]];
 805     while(current <= range.second)
 806     {
 807       this->index_pointers[i] = c;
 808       current += this->index_rate;
 809       i++;
 810     }
 811   }
 812 }
 813
 814 usint
 815 buildRanges(usint* distribution, pair_type* index_ranges)
 816 {
 817   if(distribution == 0 || index_ranges == 0) { return 0; }
 818
 819   usint sum = 0;
 820   for(usint c = 0; c < CHARS; c++)
 821   {
 822     if(distribution[c] == 0)
 823     {
 824       if(sum == 0) { index_ranges[c] = EMPTY_PAIR; }
 825       else         { index_ranges[c] = pair_type(sum, sum - 1); }
 826     }
 827     else
 828     {
 829       index_ranges[c].first = sum;
 830       sum += distribution[c];
 831       index_ranges[c].second = sum - 1;
 832     }
 833   }
 834
 835   return sum;
 836 }
 837
 838
 839 } // namespace CSA