dcover/bwt.hpp

   1 /*
   2  * bwt.hpp
   3  *
   4  * Copyright 2004 Juha Ka"rkka"inen <juha.karkkainen@cs.helsinki.fi>
   5  *
   6  */
   7
   8 #include "dcsuffixsort.hpp"
   9 #include "dcsamplesort.hpp"
  10 #include "stringsort.hpp"
  11 #include "difference_cover.hpp"
  12 #include "distribute.hpp"
  13 #include "kmpsplit.hpp"
  14 /*
  15 #include "doubling.hpp"
  16 */
  17
  18 #include <iterator>
  19 #include <vector>
  20 #include <set>
  21 #include <iostream>
  22 #include <iomanip>
  23 #include <cassert>
  24 /*
  25 #include <algorithm>
  26 #include <cstdlib>
  27 */
  28
  29 #ifndef DCOVER
  30 #define DCOVER 8
  31 #endif
  32
  33 #ifndef RANDSEED
  34 #define RANDSEED 46
  35 #endif
  36
  37 #define MIN_PIVOTS 3
  38 #define OVERSAMPLING_RATE 10.0
  39
  40
  41 struct do_nothing {
  42   template <typename Iterator>
  43   void operator() (Iterator, Iterator) {}
  44 };
  45
  46
  47
  48
  49 template <class StringIterator, class OutputIterator>
  50 typename std::iterator_traits<StringIterator>::difference_type
  51 compute_bwt (StringIterator str, StringIterator str_end,
  52              OutputIterator bwt, unsigned endmarkers, unsigned int dcover = DCOVER)
  53 {
  54
  55   typedef typename std::iterator_traits<StringIterator>::difference_type
  56                         difftype;   // signed integral type
  57   typedef unsigned long sizetype;   // unsigned integral type
  58   typedef difftype postype;   // type for representing vector positions
  59
  60   typedef typename std::vector<postype>::iterator pos_iterator;
  61
  62 #if DEBUG > 0
  63   std::cout << "BWT construction starts" << std::endl;
  64 #endif
  65
  66   //-----------------
  67   // setup and check
  68   //-----------------
  69
  70   sizetype length = std::distance(str, str_end);
  71
  72 #if DEBUG > 1
  73   typedef typename std::iterator_traits<StringIterator>::value_type
  74     chartype;
  75   typedef typename std::iterator_traits<OutputIterator>::value_type
  76     bwttype;
  77   std::cout << " text memory=" << length*sizeof(chartype) << "\n"
  78             << " bwt memory=" << (length+1)*sizeof(chartype) << std::endl;
  79 #endif
  80
  81   dc_sampler<postype> dcsampler(length+1, dcover);
  82
  83 #if DEBUG > 1
  84   std::cout << " dcsampler memory="
  85             << (dcsampler.packed_period() + 2*dcsampler.period())
  86                 * sizeof(unsigned)
  87             << std::endl;
  88 #endif
  89
  90   sizetype blocksize = dcsampler.samplesize();
  91   sizetype npivots = MIN_PIVOTS;
  92
  93   // make average bucket size <= blocksize/OVERSAMPLING_RATE
  94   while ((OVERSAMPLING_RATE*length)/(npivots+1.0) > blocksize) ++npivots;
  95
  96 #if DEBUG > 1
  97   std::cout << " length=" << length
  98             << " samplesize=" << dcsampler.samplesize()
  99             << " blocksize=" << blocksize
 100             << " npivots=" << npivots
 101             << std::endl;
 102 #endif
 103
 104   if (dcsampler.period() + 2*npivots > length) {
 105     // too small text for sampling
 106
 107 #if DEBUG > 0
 108     std::cout << "small text: sort all suffixes\n";
 109 #endif
 110
 111     std::vector<postype> sa(length);
 112     for (pos_iterator i=sa.begin(); i!=sa.end(); ++i) {
 113       *i = i - sa.begin();
 114     }
 115     suffix_mkqsort(str, str_end, sa.begin(), sa.end(), length, do_nothing());
 116     difftype eof_position = 0;
 117     for (pos_iterator i=sa.begin(); i!=sa.end(); ++i) {
 118       if (*i != 0) {
 119           if (str[*i-1] <= endmarkers)
 120               *bwt++ = 0;
 121           else
 122               *bwt++ = str[*i-1] - endmarkers;
 123       } else {
 124           if (str[*(i-1)-1] <= endmarkers)
 125               *bwt++ = 0;
 126           else
 127               *bwt++ = str[*(i-1)-1] - endmarkers;
 128           eof_position = i - sa.begin() +1; // FIXME
 129       }
 130     }
 131     return eof_position;
 132   }
 133
 134   // don't take too short suffixes as pivots
 135   unsigned int pivotrange = length-dcsampler.period()+1;
 136
 137   //-----------------------------------
 138   // setup difference cover sample
 139   //-----------------------------------
 140
 141 #if DEBUG > 0
 142   std::cout << "difference cover construction ... " << std::endl;
 143 #endif
 144
 145   sizetype samplesize = dcsampler.samplesize();
 146   std::vector<difftype> dcranks(samplesize);
 147   std::vector<difftype> suffixes(samplesize);
 148
 149 #if DEBUG > 1
 150   std::cout << "dcranks memory="
 151             << samplesize * sizeof(difftype)
 152             << "\n";
 153   std::cout << "suffixes memory="
 154             << samplesize * sizeof(difftype)
 155             << std::endl;
 156 #endif
 157
 158   sort_dc_sample(dcsampler, str, str_end,
 159                  dcranks.begin(), dcranks.end(),
 160                  suffixes.begin(), suffixes.end());
 161
 162 #if DEBUG > 3
 163   std::cout << "sample ranks:\n";
 164   std::copy(dcranks.begin(), dcranks.end(),
 165             std::ostream_iterator<int>(std::cout,"\n"));
 166   std::cout << "end of sample ranks" << std::endl;
 167 #endif
 168
 169
 170
 171   //-------------------
 172   // setup pivots
 173   //-------------------
 174
 175 #if DEBUG > 0
 176   std::cout << "setting up pivots ..." << std::endl;
 177 #endif
 178
 179   std::set<difftype> pivotset;
 180   double scale_factor = (pivotrange * 1.0L)/(RAND_MAX+1.0L);
 181   std::srand(RANDSEED);
 182   while (pivotset.size() < npivots) {
 183     int r = std::rand();
 184     int newpivot = static_cast<difftype>(scale_factor*r);
 185     pivotset.insert(newpivot);
 186   }
 187
 188 #if DEBUG > 1
 189   std::cout << " pivotset memory = " << pivotset.size()
 190             << " elements" << std::endl;
 191 #endif
 192
 193   std::vector<difftype> pivots(pivotset.begin(), pivotset.end());
 194   pivotset.clear();
 195
 196   // sort the pivots
 197   dc_sort_suffixes(str, str_end, pivots.begin(), pivots.end(),
 198                    dcsampler, dcranks.begin(), dcranks.end());
 199
 200 #if DEBUG > 2
 201   std::cout << "pivots sorted lexicographically:\n";
 202   for (pos_iterator i = pivots.begin(); i != pivots.end(); ++i) {
 203     std::cout << std::setw(8) << *i << " ";
 204     std::copy(str+*i, std::min(str_end, str+*i+40),
 205               std::ostream_iterator<char>(std::cout,""));
 206     std::cout << std::endl;
 207   }
 208 #endif
 209
 210
 211 #if DEBUG > 0
 212   std::cout << " constructing distributor ... " << std::endl;
 213 #endif
 214
 215   suffix_distributor<StringIterator, pos_iterator,
 216                      dc_sampler<postype>, pos_iterator>
 217     distributor(str, str_end, pivots.begin(), pivots.end(),
 218                 dcsampler, dcranks.begin());
 219
 220 #if DEBUG > 1
 221   distributor.print();
 222 #endif
 223
 224
 225
 226   //------------------------
 227   // compute bucketsizes
 228   //------------------------
 229
 230 #if DEBUG > 0
 231   std::cout << "computing bucket sizes ... " << std::endl;
 232 #endif
 233
 234   std::vector<difftype> buckets(distributor.nbuckets());
 235
 236 #if DEBUG > 1
 237   std::cout << " number of buckets: " << distributor.nbuckets() << "\n";
 238   std::cout << " bucket vector memory="
 239             << buckets.size() * sizeof(difftype)
 240             << std::endl;
 241 #endif
 242
 243   for (StringIterator suf = str; suf != str_end; ++suf) {
 244     unsigned int bucket = distributor.find_bucket(suf);
 245     ++buckets[bucket];
 246
 247 #if DEBUG > 3
 248     std::cout << " into bucket " << bucket;
 249     std::cout << ":" << std::setw(7) << suf-str << " ";
 250     std::copy(suf, std::min(str_end, suf+40),
 251               std::ostream_iterator<char>(std::cout,""));
 252     std::cout << std::endl;
 253 #endif
 254
 255   }
 256   // ++buckets[0]; // empty suffix
 257
 258 #if DEBUG > 1
 259   std::cout << " bucket sizes: ";
 260   std::copy(buckets.begin(), buckets.end(),
 261             std::ostream_iterator<int>(std::cout," "));
 262   std::cout << std::endl;
 263 #endif
 264
 265
 266   //-----------------------
 267   // build and sort blocks
 268   //-----------------------
 269
 270 #if DEBUG > 0
 271   std::cout << "Computing and sorting blocks ... " << std::endl;
 272 #endif
 273
 274 //  *bwt++ = *(str_end-1);  //FIXME no empty suffix?
 275
 276   postype eof_position = 0;
 277   sizetype total_sum = 1;
 278
 279   postype bucket_begin = 0;
 280   while (bucket_begin < (postype)buckets.size()) {
 281
 282     sizetype sum = buckets[bucket_begin];
 283     postype bucket_end = bucket_begin + 1;
 284     while (bucket_end < (postype)buckets.size()
 285            && sum+buckets[bucket_end] <= blocksize) {
 286       sum += buckets[bucket_end];
 287       ++bucket_end;
 288     }
 289
 290     postype left_pivot = bucket_begin-1;
 291     postype right_pivot = bucket_end-1;
 292
 293     short lcp = 0;
 294     if (left_pivot >= 0 && right_pivot < (postype)pivots.size()) {
 295 #if __WORDSIZE == 32
 296       lcp = suffix_lcp(str+pivots[left_pivot], str+pivots[right_pivot],
 297                        str_end, 0l, dcsampler.period()-1);
 298 #else
 299       lcp = suffix_lcp(str+pivots[left_pivot], str+pivots[right_pivot],
 300                str_end, 0l, dcsampler.period()-1);
 301 #endif
 302     }
 303
 304 #if DEBUG > 1
 305     std::cout << "block [" << bucket_begin << "," << bucket_end << ")"
 306               << " size=" << sum
 307               << " lcp=" << lcp
 308               << std::endl;
 309 #endif
 310
 311     if (sum > blocksize) {
 312       blocksize = sum;
 313       suffixes.resize(blocksize);
 314
 315 #if DEBUG > 1
 316       std::cout << " blocksize increased to "
 317                 << blocksize << "\n";
 318       std::cout << " suffixes memory resized to "
 319                 << blocksize*sizeof(postype)
 320                 << std::endl;
 321 #endif
 322
 323     }
 324
 325     pos_iterator i = suffixes.begin();
 326     //if (bucket_begin == 0) { *i++ = length;       // empty suffix }
 327
 328     postype left_pivot_not_beyond = left_pivot<0 ? 0 : left_pivot;
 329     kmp_split<StringIterator, dc_sampler<postype>, pos_iterator>
 330       left_cmp(str, str_end, str+pivots[left_pivot_not_beyond],
 331                    dcsampler, dcranks.begin());
 332
 333     postype right_pivot_not_beyond =
 334       right_pivot==(postype)pivots.size() ? right_pivot-1 : right_pivot;
 335     kmp_split<StringIterator, dc_sampler<postype>, pos_iterator>
 336       right_cmp(str, str_end, str+pivots[right_pivot_not_beyond],
 337                    dcsampler, dcranks.begin());
 338
 339     for (StringIterator suf = str; suf != str_end; ++suf) {
 340
 341 #if DEBUG > 3
 342       std::cout << ":" << std::setw(3) << suf-str << " ";
 343       std::copy(suf, std::min(str_end, suf+40),
 344                 std::ostream_iterator<char>(std::cout,""));
 345       std::cout << std::endl;
 346 #endif
 347
 348       /*
 349       // is the suffix left of the block
 350       if (left_pivot >=0) {
 351         postype pivotpos = pivots[left_pivot];
 352         if (*suf < str[pivotpos]) continue;
 353         int lcp = suffix_lcp(suf, str+pivotpos,
 354                              str_end, 0, dcsampler.period()-1);
 355         if (lcp < dcsampler.period()-1) {
 356           if (suf+lcp == str_end || suf[lcp] < (str+pivotpos)[lcp]) {
 357             continue;
 358           }
 359         } else {
 360           postype sufpos = suf-str;
 361           std::pair<postype,postype> result
 362             = dcsampler.pair(sufpos, pivotpos);
 363           if (dcranks[result.first] < dcranks[result.second]){
 364             continue;
 365           }
 366         }
 367       }
 368
 369       // is the suffix right of the block
 370       if (right_pivot < (postype)pivots.size()) {
 371         postype pivotpos = pivots[right_pivot];
 372         if (*suf > str[pivotpos]) continue;
 373         int lcp = suffix_lcp(str+pivotpos, suf,
 374                              str_end, 0, dcsampler.period()-1);
 375         if (lcp < dcsampler.period()-1) {
 376           if (suf+lcp < str_end && (str+pivotpos)[lcp] < suf[lcp]) {
 377             continue;
 378           }
 379         } else {
 380           postype sufpos = suf-str;
 381           std::pair<postype,postype> result
 382             = dcsampler.pair(pivotpos, sufpos);
 383           if (dcranks[result.first] <= dcranks[result.second]){
 384             continue;
 385           }
 386         }
 387       }
 388       */
 389
 390       //postype bucket = distributor.find_bucket(suf);
 391
 392       bool to_left = left_pivot>0 && left_cmp.is_next_less();
 393       bool to_right = right_pivot<(postype)pivots.size() && !right_cmp.is_next_less() ;
 394
 395       if (to_left || to_right) {
 396         //assert (bucket_begin > bucket || bucket >= bucket_end);
 397         continue;
 398       }
 399
 400       // it is in the bucket
 401       //postype bucket = distributor.find_bucket(suf);
 402
 403       //if (bucket_begin > bucket || bucket >= bucket_end) continue;
 404
 405       //assert (bucket_begin <= bucket && bucket < bucket_end);
 406
 407       *i++ = suf-str;
 408 #if DEBUG > 3
 409       //std::cout << "into bucket " << bucket;
 410       std::cout << ":" << std::setw(3) << suf-str << " ";
 411       std::copy(suf, std::min(str_end, suf+40),
 412                 std::ostream_iterator<char>(std::cout,""));
 413       std::cout << std::endl;
 414 #endif
 415
 416     }
 417
 418     assert(i == suffixes.begin()+sum);
 419
 420     dc_sort_suffixes(str, str_end, suffixes.begin(), i,
 421                      dcsampler, dcranks.begin(), dcranks.end());
 422
 423 #if DEBUG > 3
 424     std::cout << "sorted bucket:\n";
 425     for (pos_iterator j = suffixes.begin(); j != i; ++j) {
 426       std::cout << std::setw(3) << *j << " ";
 427       std::copy(str+*j, std::min(str_end, str+*j+40),
 428                 std::ostream_iterator<char>(std::cout,""));
 429       std::cout << std::endl;
 430     }
 431 #endif
 432
 433     for (pos_iterator j=suffixes.begin(); j != i; ++j) {
 434       if (*j != 0) {
 435           if (str[*j-1] <= endmarkers)
 436               *bwt++ = 0;
 437           else
 438               *bwt++ = str[*j-1] - endmarkers;
 439       } else {
 440           if (str[*(j-1)-1] <= endmarkers)
 441               *bwt++ = 0;
 442           else
 443               *bwt++ = str[*(j-1)-1] - endmarkers;
 444         eof_position = total_sum + (j - suffixes.begin());
 445       }
 446     }
 447
 448     total_sum += sum;
 449
 450     bucket_begin = bucket_end;
 451   }
 452
 453   assert(total_sum == length+1);
 454
 455   return eof_position;
 456 }