incbwt/rlcsa.h

   1 #ifndef RLCSA_H
   2 #define RLCSA_H
   3
   4 #include <fstream>
   5 #include <vector>
   6
   7 #include "bits/deltavector.h"
   8 #include "bits/rlevector.h"
   9 #include "bits/nibblevector.h"
  10
  11 #include "sasamples.h"
  12 #include "lcpsamples.h"
  13 #include "misc/parameters.h"
  14
  15
  16 namespace CSA
  17 {
  18
  19
  20 const std::string SA_EXTENSION = ".sa";
  21 const std::string PSI_EXTENSION = ".psi";
  22 const std::string ARRAY_EXTENSION = ".rlcsa.array";
  23 const std::string SA_SAMPLES_EXTENSION = ".rlcsa.sa_samples";
  24 const std::string PARAMETERS_EXTENSION = ".rlcsa.parameters";
  25 const std::string LCP_SAMPLES_EXTENSION = ".lcp_samples";
  26 const std::string PLCP_EXTENSION = ".plcp";
  27
  28
  29 const parameter_type RLCSA_BLOCK_SIZE  = parameter_type("RLCSA_BLOCK_SIZE", 32);
  30 const parameter_type SAMPLE_RATE       = parameter_type("SAMPLE_RATE", 512);
  31 const parameter_type SUPPORT_LOCATE    = parameter_type("SUPPORT_LOCATE", 1);
  32 const parameter_type SUPPORT_DISPLAY   = parameter_type("SUPPORT_DISPLAY", 1);
  33
  34
  35 #ifdef USE_NIBBLE_VECTORS
  36 typedef NibbleVector  PsiVector;
  37 #else
  38 typedef RLEVector     PsiVector;
  39 #endif
  40
  41
  42 class RLCSA
  43 {
  44   public:
  45
  46 //--------------------------------------------------------------------------
  47 //  CONSTRUCTION
  48 //--------------------------------------------------------------------------
  49
  50     const static usint ENDPOINT_BLOCK_SIZE = 16;
  51
  52     RLCSA(const std::string& base_name, bool print = false);
  53
  54     /*
  55       If multiple_sequences is true, each 0 is treated as a end of sequence marker.
  56       There must be nonzero characters between the 0s. data[bytes - 1] must also be 0.
  57       FIXME Crashes if bytes >= 2 GB.
  58     */
  59     RLCSA(uchar* data, usint bytes, usint block_size, usint sa_sample_rate, bool multiple_sequences, bool delete_data);
  60
  61     // Destroys contents of index and increment.
  62     // threads has no effect unless MULTITHREAD_SUPPORT is defined
  63     RLCSA(RLCSA& index, RLCSA& increment, usint* positions, usint block_size, usint threads = 1);
  64     ~RLCSA();
  65
  66     void writeTo(const std::string& base_name) const;
  67
  68     bool isOk() const;
  69
  70 //--------------------------------------------------------------------------
  71 //  QUERIES
  72 //--------------------------------------------------------------------------
  73
  74     // Returns the closed interval of suffix array containing the matches.
  75     pair_type count(const std::string& pattern) const;
  76
  77     // Used when merging CSAs.
  78     void reportPositions(uchar* data, usint length, usint* positions) const;
  79
  80     // Returns SA[range]. User must free the buffer. Latter version uses buffer provided by the user.
  81     usint* locate(pair_type range) const;
  82     usint* locate(pair_type range, usint* data) const;
  83
  84     // Returns SA[index].
  85     usint locate(usint index) const;
  86
  87     // Returns T^{sequence}[range]. User must free the buffer.
  88     // Third version uses buffer provided by the user.
  89     uchar* display(usint sequence) const;
  90     uchar* display(usint sequence, pair_type range) const;
  91     uchar* display(usint sequence, pair_type range, uchar* data) const;
  92
  93     // Displays the intersection of T[position - context, position + len + context - 1]
  94     // and T^{getSequenceForPosition(position)}.
  95     // This is intended for displaying an occurrence of a pattern of length 'len'
  96     // with 'context' extra characters on both sides.
  97     // The actual length of the returned string is written into result_length.
  98     uchar* display(usint position, usint len, usint context, usint& result_length) const;
  99
 100     // Returns the actual length of the prefix. User must provide the buffer.
 101     usint displayPrefix(usint sequence, usint len, uchar* data) const;
 102
 103     // Returns at most max_len characters starting from T[SA[index]].
 104     // User must provide the buffer. Returns the number of characters in buffer.
 105     usint displayFromPosition(usint index, usint max_len, uchar* data) const;
 106
 107     // Get the range of SA values for the sequence identified by
 108     // a sequence number or a SA value.
 109     pair_type getSequenceRange(usint number) const;
 110     pair_type getSequenceRangeForPosition(usint value) const;
 111
 112     // Get the sequence number for given SA value(s).
 113     // The returned array is the same as the parameter.
 114     usint getSequenceForPosition(usint value) const;
 115     usint* getSequenceForPosition(usint* value, usint length) const;
 116
 117     // Changes SA value to (sequence, offset).
 118     pair_type getRelativePosition(usint value) const;
 119
 120     // Returns the BWT of the collection including end of sequence markers.
 121     uchar* readBWT() const;
 122     uchar* readBWT(pair_type range) const;
 123
 124     // Returns the number of equal letter runs in the BWT. Runs consisting of end markers are ignored.
 125     usint countRuns() const;
 126
 127 //--------------------------------------------------------------------------
 128 //  BASIC OPERATIONS
 129 //--------------------------------------------------------------------------
 130
 131     // The return values of these functions include the implicit end markers.
 132     // To get SA indexes, subtract getNumberOfSequences() from the value.
 133
 134     inline usint psi(usint index) const
 135     {
 136       if(index >= this->data_size)
 137       {
 138         return this->data_size + this->number_of_sequences;
 139       }
 140
 141       usint c = this->getCharacter(index);
 142       return this->psiUnsafe(index, c);
 143     }
 144
 145     // This version returns a run.
 146     inline pair_type psi(usint index, usint max_length) const
 147     {
 148       if(index >= this->data_size)
 149       {
 150         return pair_type(this->data_size + this->number_of_sequences, 0);
 151       }
 152
 153       usint c = this->getCharacter(index);
 154       PsiVector::Iterator iter(*(this->array[c]));
 155       return iter.selectRun(index - this->index_ranges[c].first, max_length);
 156     }
 157
 158     inline usint C(usint c) const
 159     {
 160       if(c >= CHARS)
 161           return this->data_size + this->number_of_sequences;
 162       return this->index_ranges[c].first + this->number_of_sequences - 1;
 163     }
 164
 165     inline usint LF(usint index, usint c) const
 166     {
 167       if(c >= CHARS)
 168       {
 169         return this->data_size + this->number_of_sequences;
 170       }
 171       if(this->array[c] == 0)
 172       {
 173         if(c < this->text_chars[0]) { return this->number_of_sequences - 1; }
 174         return this->index_ranges[c].first + this->number_of_sequences - 1;
 175       }
 176       index += this->number_of_sequences;
 177
 178       PsiVector::Iterator iter(*(this->array[c]));
 179       return this->LF(index, c, iter);
 180     }
 181
 182 //--------------------------------------------------------------------------
 183 //  REPORTING
 184 //--------------------------------------------------------------------------
 185
 186     inline bool supportsLocate() const { return this->support_locate; }
 187     inline bool supportsDisplay() const { return this->support_display; }
 188     inline usint getSize() const { return this->data_size; }
 189     inline usint getNumberOfSequences() const { return this->number_of_sequences; }
 190     inline usint getBlockSize() const { return this->array[this->text_chars[0]]->getBlockSize(); }
 191
 192     // Returns the size of the data structure.
 193     usint reportSize(bool print = false) const;
 194
 195     void printInfo() const;
 196
 197 //--------------------------------------------------------------------------
 198 //  LCP EXPERIMENTS
 199 //--------------------------------------------------------------------------
 200
 201     // Optimized version:
 202     //   - Interleaves main loop with computing irreducible values.
 203     //   - Encodes maximal runs from a true local maximum to a true local minimum.
 204     RLEVector* buildPLCP(usint block_size) const;
 205
 206     // Returns the number of samples. sampled_values will be a pointer to the samples.
 207     usint sampleLCP(usint sample_rate, pair_type*& sampled_values, bool report = false) const;
 208
 209     usint lcp(usint index, const LCPSamples& lcp_samples) const;
 210     usint lcp(usint index, const LCPSamples& lcp_samples, const RLEVector& plcp) const;
 211
 212     // Calculate LCP[index] directly.
 213     usint lcpDirect(usint index) const;
 214
 215     // Writes PLCP[start] to PLCP[stop - 1].
 216     inline void encodePLCPRun(RLEEncoder& plcp, usint start, usint stop, usint first_val) const
 217     {
 218       plcp.setRun(2 * start + first_val, stop - start);
 219 //      std::cerr << "(" << start << ", " << stop << ", " << first_val << ")" << std::endl;
 220     }
 221
 222 //--------------------------------------------------------------------------
 223 //  INTERNAL VARIABLES
 224 //--------------------------------------------------------------------------
 225
 226   private:
 227     bool ok;
 228
 229     PsiVector* array[CHARS];
 230     SASamples* sa_samples;
 231
 232     pair_type index_ranges[CHARS];
 233     usint data_size;
 234
 235     usint text_chars[CHARS];  // which characters are present in the text
 236     usint chars;
 237
 238     usint index_pointers[CHARS]; // which of the above is at i * index_rate
 239     usint index_rate;
 240
 241     bool support_locate, support_display;
 242     usint sample_rate;
 243
 244     usint number_of_sequences;
 245     DeltaVector* end_points;
 246
 247 //--------------------------------------------------------------------------
 248 //  INTERNAL VERSIONS OF QUERIES
 249 //--------------------------------------------------------------------------
 250
 251     void locateUnsafe(pair_type range, usint* data) const;
 252     bool processRun(pair_type run, usint* data, usint* offsets, bool* finished, PsiVector::Iterator** iters) const;
 253     void displayUnsafe(pair_type range, uchar* data) const;
 254
 255 //--------------------------------------------------------------------------
 256 //  INTERNAL VERSIONS OF BASIC OPERATIONS
 257 //--------------------------------------------------------------------------
 258
 259     inline usint psi(usint index, PsiVector::Iterator** iters) const
 260     {
 261       usint c = this->getCharacter(index);
 262       return iters[c]->select(index - this->index_ranges[c].first);
 263     }
 264
 265     inline pair_type psi(usint index, usint max_length, PsiVector::Iterator** iters) const
 266     {
 267       usint c = this->getCharacter(index);
 268       return iters[c]->selectRun(index - this->index_ranges[c].first, max_length);
 269     }
 270
 271     // Returns psi(index), assuming the suffix of rank index begins with c.
 272     inline usint psiUnsafe(usint index, usint c) const
 273     {
 274       PsiVector::Iterator iter(*(this->array[c]));
 275       return this->psiUnsafe(index, c, iter);
 276     }
 277
 278     // As above, but with a given iterator.
 279     inline usint psiUnsafe(usint index, usint c, PsiVector::Iterator& iter) const
 280     {
 281       return iter.select(index - this->index_ranges[c].first);
 282     }
 283
 284     // As above, but returns the next value of psi.
 285     inline usint psiUnsafeNext(usint c, PsiVector::Iterator& iter) const
 286     {
 287       return iter.selectNext();
 288     }
 289
 290     inline pair_type backwardSearchStep(pair_type range, usint c) const
 291     {
 292       if(this->array[c] == 0) { return EMPTY_PAIR; }
 293       PsiVector::Iterator iter(*(this->array[c]));
 294
 295       usint start = this->index_ranges[c].first + this->number_of_sequences - 1;
 296       range.first = start + iter.rank(range.first, true);
 297       range.second = start + iter.rank(range.second);
 298
 299       return range;
 300     }
 301
 302     inline usint LF(usint index, usint c, PsiVector::Iterator& iter) const
 303     {
 304       return this->index_ranges[c].first + this->number_of_sequences + iter.rank(index) - 1;
 305     }
 306
 307 //--------------------------------------------------------------------------
 308 //  INTERNAL STUFF
 309 //--------------------------------------------------------------------------
 310
 311     // Creates an array of iterators for every vector in this->array.
 312     PsiVector::Iterator** getIterators() const;
 313     void deleteIterators(PsiVector::Iterator** iters) const;
 314
 315     inline usint getCharacter(usint pos) const
 316     {
 317       const usint* curr = &(this->text_chars[this->index_pointers[pos / this->index_rate]]);
 318       while(pos > this->index_ranges[*curr].second) { curr++; }
 319       return *curr;
 320     }
 321
 322     void mergeEndPoints(RLCSA& index, RLCSA& increment);
 323     void mergeSamples(RLCSA& index, RLCSA& increment, usint* positions);
 324
 325     void buildCharIndexes(usint* distribution);
 326
 327     // These are not allowed.
 328     RLCSA();
 329     RLCSA(const RLCSA&);
 330     RLCSA& operator = (const RLCSA&);
 331 };
 332
 333
 334 // Returns the total number of characters.
 335 usint buildRanges(usint* distribution, pair_type* index_ranges);
 336
 337
 338 } // namespace CSA
 339
 340
 341 #endif // RLCSA_H