RLWaveletTree.cpp

   1 /***************************************************************************
   2  *   Copyright (C) 2007 by Veli Mäkinen                                    *
   3  *                                                                         *
   4  *                                                                         *
   5  *   This program is free software; you can redistribute it and/or modify  *
   6  *   it under the terms of the GNU General Public License as published by  *
   7  *   the Free Software Foundation; either version 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  *                                                                         *
  10  *   This program is distributed in the hope that it will be useful,       *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  13  *   GNU General Public License for more details.                          *
  14  *                                                                         *
  15  *   You should have received a copy of the GNU General Public License     *
  16  *   along with this program; if not, write to the                         *
  17  *   Free Software Foundation, Inc.,                                       *
  18  *   51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.         *
  19  ***************************************************************************/
  20
  21 #include "RLWaveletTree.h"
  22 #include <stdexcept>
  23 using namespace std;
  24
  25 void RLWaveletTree::remapAlphabet(uchar *text, ulong n)
  26 {
  27     ulong *chars = new ulong[256/W+1];
  28     ulong i;
  29     for (i = 0; i < 256/W+1; ++i)
  30         chars[i] = 0;
  31     // Set up flags:
  32     for (i = 0; i < n; ++i)
  33         Tools::SetField(chars, 1, text[i], 1);
  34     charMap = new BitRank(chars, 256, true);
  35
  36     // Remap text
  37     for (i = 0; i < n; ++i)
  38         text[i] = charMap->rank(text[i])-1;
  39 }
  40
  41 RLWaveletTree::RLWaveletTree(uchar *D, ulong n) {
  42     remapAlphabet(D, n);
  43 //    std::cout << "Number of chars: " << charMap->rank(255) << ", maxlevel = " << Tools::CeilLog2(charMap->rank(255)) << std::endl;
  44
  45     this->n = n;
  46     this->maxlevel = Tools::CeilLog2(charMap->rank(255));
  47
  48     ulong i, k, value;
  49     ulong** bucket = new ulong*[maxlevel+1];
  50     ulong *wtree = new ulong[n*maxlevel/(8*sizeof(ulong))+1];
  51     for (i = 0; i < n*maxlevel/(8*sizeof(ulong))+1; i ++)
  52         wtree[i] = 0lu;
  53
  54    // allocate memory for storing the locations of boudaries
  55     for (k=0; k<=maxlevel; k++) {
  56         bucket[k] = (ulong *)new ulong[(ulong)1 << k];
  57         for (i=0; i<((ulong)1 << k); i++)
  58             bucket[k][i]=0;
  59     }
  60
  61    // compute the sizes of buckets
  62     for (i=0; i<n; i++) {
  63         value = D[i];
  64         for (k=0; k<=maxlevel; k++)
  65             bucket[k][value >> (maxlevel-k)]++;
  66     }
  67
  68    // cumulate the bucket sizes into boundary locations
  69     for (k=0; k<=maxlevel; k++) {
  70         for (i=1; i<((ulong)1 << k); i++)
  71             bucket[k][i]+=bucket[k][i-1];
  72         for (i=((ulong)1 << k)-1; i>=1; i--)
  73             bucket[k][i]=bucket[k][i-1];
  74         bucket[k][0]=0;
  75     }
  76
  77    //store the leaf information
  78     leaves = new ulong[((ulong)1 << maxlevel)+1];
  79     for (k=0; k<(ulong)1<<maxlevel; k++)
  80         leaves[k] = bucket[maxlevel][k];
  81     leaves[k]=n;
  82
  83    // compute wtree using the information about boundaries
  84     for (i=0; i<n; i++) {
  85         value = D[i];
  86         for (k=0; k<maxlevel; k++)
  87             Tools::SetField(wtree,1,k*n+(bucket[k][value >> (maxlevel-k)])++,(value>>(maxlevel-k-1)) & (ulong)1);
  88     }
  89     delete [] D;
  90     // deallocate boundaries
  91     for (k=0; k<=maxlevel; k++)
  92         delete[] bucket[k];
  93     delete[] bucket;
  94
  95     bitrank = new GapEncode(wtree, n*maxlevel, true);
  96     prevRunLen = 0;
  97     prevChar = 0;
  98     prevPos = 0;
  99 }
 100
 101 RLWaveletTree::~RLWaveletTree() {
 102     delete charMap;
 103     delete [] leaves;
 104     delete bitrank;
 105 }
 106
 107 ulong RLWaveletTree::rank(ulong c, ulong i) {
 108 //    if (i >= n) FIXME Replace with assert?
 109 //        std::cout << "RLWT::rank parameter i out of range" << std::endl;
 110    ulong k;
 111    ulong j=i;
 112    ulong left = 0, right = (ulong)1 << maxlevel;
 113    ulong before=0;
 114    c = charMap->rank(c)-1;
 115
 116    // first round separately to avoid extra "if"
 117    if ((c >> (maxlevel-1)) & (ulong)1) { // go to the right
 118       j = bitrank->rank(j)-1;
 119       left = right>>1;
 120    }
 121    else { // go to the left
 122       j -= bitrank->rank(j);
 123       right = right>>1;
 124    }
 125    // then the next without the extra if
 126    for (k=1;k<maxlevel;k++) {
 127       before = bitrank->rank(n*k+leaves[left]-1);
 128       if ((c >> (maxlevel-k-1)) & (ulong)1) { // go to the right
 129          j = bitrank->rank(n*k+leaves[left]+j)-before-1;
 130          left = (right+left)>>1;
 131       }
 132       else { // go to the left
 133          j -= bitrank->rank(n*k+leaves[left]+j)-before;
 134          right = (right+left)>>1;
 135       }
 136       // can never happen: if (left > right) return 0;
 137    }
 138    return j+1;
 139 }
 140
 141 ulong RLWaveletTree::select(ulong c, ulong j) {
 142     c = charMap->rank(c)-1;
 143    ulong left=c, right=c+1;
 144    ulong k, i=j-1, before;
 145    for (k=maxlevel-1; k>0;--k) {
 146        if ((c >> (maxlevel-k-1)) & (ulong)1) { // right side of parent
 147           left -= right-left;
 148           before = bitrank->rank(n*k+leaves[left]-1);
 149           i = bitrank->select(before+i+1)-n*k-leaves[left];
 150        }
 151        else { // left side of parent
 152           right +=right-left;
 153           before = n*k+leaves[left]-bitrank->rank(n*k+leaves[left]-1);
 154           i = bitrank->select0(before+i+1)-n*k-leaves[left];
 155        }
 156    }
 157    // last level separately to avoid one extra if
 158    if ((c >> (maxlevel-1)) & (ulong)1) // right side of parent
 159       return bitrank->select(i+1);
 160    else // left side of parent
 161       return bitrank->select0(i+1);
 162 }
 163
 164 bool RLWaveletTree::IsCharAtPos(ulong c, ulong i) {
 165    return charAtPos(i)==c;
 166 }
 167
 168
 169 ulong RLWaveletTree::charAtPos(ulong i)
 170 {
 171 //    if (i >= n) // FIXME Replace with assert?
 172 //        std::cout << "RLWT::charAtPos parameter i out of range" << std::endl;
 173     ulong k;
 174     ulong j=i;
 175     ulong left = 0, right = (ulong)1 << maxlevel;
 176     ulong before=0;
 177     ulong c = 0;
 178     ulong rank_tmp = 0;
 179     // first round separately to avoid extra "if"
 180     if (bitrank->IsBitSet(j, &rank_tmp)) { // go to the right
 181         c = c | ((ulong)1 << (maxlevel-1));
 182         j = rank_tmp-1; //bitrank->rank(j)-1;
 183         left = right>>1;
 184     }
 185     else { // go to the left
 186         // c = c
 187         j -= rank_tmp; // bitrank->rank(j);
 188         right = right>>1;
 189     }
 190     // then the next without the extra if
 191     for (k=1;k<maxlevel;k++) {
 192         before = bitrank->rank(n*k+leaves[left]-1);
 193         if (bitrank->IsBitSet(n*k+leaves[left]+j, &rank_tmp)) { // go to the right
 194             c = c | ((ulong)1 << (maxlevel-k-1));
 195             j = rank_tmp-before-1; //bitrank->rank(n*k+leaves[left]+j)-before-1;
 196             left = (right+left)>>1;
 197         }
 198         else { // go to the left
 199             // c = c
 200             j -= rank_tmp-before; //bitrank->rank(n*k+leaves[left]+j)-before;
 201             right = (right+left)>>1;
 202         }
 203         // can never happen: if (left > right) return 0;
 204     }
 205     return charMap->select(c+1);
 206 }
 207
 208 /* returns also the rank-1 of char at given position i */
 209 ulong RLWaveletTree::charAtPos(ulong i, ulong *rank)
 210 {
 211 //    if (i >= n) // FIXME Replace with assert?
 212 //        std::cout << "RLWT::charAtPos2 parameter i out of range" << std::endl;
 213     ulong k;
 214     ulong j=i;
 215     ulong left = 0, right = (ulong)1 << maxlevel;
 216     ulong before=0;
 217     ulong c = 0;
 218     ulong rank_tmp = 0;
 219     // first round separately to avoid extra "if"
 220     if (bitrank->IsBitSet(j, &rank_tmp)) { // FIXME j or leaves[left]+j ?
 221         c = c | ((ulong)1 << (maxlevel-1));
 222         j = rank_tmp-1; //bitrank->rank(j)-1;
 223         left = right>>1;
 224     }
 225     else { // go to the left
 226         // c = c
 227         j -= rank_tmp; //bitrank->rank(j);
 228         right = right>>1;
 229     }
 230     // then the next without the extra if
 231     for (k=1;k<maxlevel;k++) {
 232         before = bitrank->rank(n*k+leaves[left]-1);
 233         if (bitrank->IsBitSet(n*k+leaves[left]+j, &rank_tmp)) { // go to the right
 234             c = c | ((ulong)1 << (maxlevel-k-1));
 235             j = rank_tmp-before-1; //bitrank->rank(n*k+leaves[left]+j)-before-1;
 236             left = (right+left)>>1;
 237         }
 238         else { // go to the left
 239             // c = c
 240             j -= rank_tmp-before; //bitrank->rank(n*k+leaves[left]+j)-before;
 241             right = (right+left)>>1;
 242         }
 243         // can never happen: if (left > right) return 0;
 244     }
 245     *rank = j;
 246     return charMap->select(c+1);
 247 }
 248
 249 /* returns also the rank-1 of char at given position i
 250  * and trailing char-run length */
 251 ulong RLWaveletTree::charAtPosRun(ulong i, ulong *rank)
 252 {
 253     // Check for buffered data from previous call
 254     if (i > prevPos && i <= prevPos + prevRunLen)
 255     {
 256         *rank = prevRank + (i-prevPos);
 257         return prevChar;
 258     }
 259
 260     prevPos = i;
 261 //    if (i >= n) // FIXME Replace with assert?
 262 //        std::cout << "RLWT::charAtPos2 parameter i out of range" << std::endl;
 263     ulong k;
 264     ulong j=i;
 265     ulong left = 0, right = (ulong)1 << maxlevel;
 266     ulong before=0;
 267     ulong c = 0;
 268     ulong rank_tmp = 0;
 269     ulong runl_tmp = 0;
 270     // first round separately to avoid extra "if"
 271     if (bitrank->IsBitSet(j, &rank_tmp, &prevRunLen)) {
 272         if (j+prevRunLen > n)
 273             prevRunLen = n-j-1;
 274         c = c | ((ulong)1 << (maxlevel-1));
 275         j = rank_tmp-1; //bitrank->rank(j)-1;
 276         left = right>>1;
 277     }
 278     else { // go to the left
 279         // c = c
 280         if (j+prevRunLen > n)
 281             prevRunLen = n-j-1;
 282         j -= rank_tmp; //bitrank->rank(j);
 283         right = right>>1;
 284     }
 285
 286     // then the next without the extra if
 287     for (k=1;k<maxlevel;k++) {
 288         before = bitrank->rank(n*k+leaves[left]-1);
 289         if (bitrank->IsBitSet(n*k+leaves[left]+j, &rank_tmp, &runl_tmp)) { // go to the right
 290             if (leaves[left]+j+runl_tmp > leaves[right])
 291                 runl_tmp = leaves[right] - (leaves[left]+j) - 1;
 292             if (prevRunLen > runl_tmp) // FIXME Better way to keep track of a run?
 293                 prevRunLen = runl_tmp;
 294             c = c | ((ulong)1 << (maxlevel-k-1));
 295             j = rank_tmp-before-1; //bitrank->rank(n*k+leaves[left]+j)-before-1;
 296             left = (right+left)>>1;
 297         }
 298         else { // go to the left
 299             if (leaves[left]+j+runl_tmp > leaves[right])
 300                 runl_tmp = leaves[right] - (leaves[left]+j) - 1;
 301             if (prevRunLen > runl_tmp)  // FIXME Better way to keep track of a run?
 302                 prevRunLen = runl_tmp;
 303             // c = c
 304             j -= rank_tmp-before; //bitrank->rank(n*k+leaves[left]+j)-before;
 305             right = (right+left)>>1;
 306         }
 307         // can never happen: if (left > right) return 0;
 308     }
 309     *rank = j;
 310     prevRank = j;
 311     prevChar = charMap->select(c+1);
 312     return charMap->select(c+1);
 313 }
 314
 315 /**
 316  * Saving data fields:
 317         ulong n;
 318         ulong maxlevel;
 319         ulong *leaves; // stores the leaf positions of characters 0..2^maxlevel-1
 320         GapEncode *bitrank;
 321         BitRank *charMap;
 322 */
 323 void RLWaveletTree::Save(FILE *file) const
 324 {
 325     if (std::fwrite(&(this->n), sizeof(ulong), 1, file) != 1)
 326         throw std::runtime_error("RLWaveletTree::Save(): file write error (n).");
 327     if (std::fwrite(&(this->maxlevel), sizeof(ulong), 1, file) != 1)
 328         throw std::runtime_error("RLWaveletTree::Save(): file write error (maxlevel).");
 329
 330     for (ulong offset = 0; offset < ((ulong)1 << maxlevel)+1; ++offset)
 331     {
 332         if (std::fwrite(this->leaves + offset, sizeof(ulong), 1, file) != 1)
 333             throw std::runtime_error("RLWaveletTree::Save(): file write error (leaves).");
 334     }
 335
 336     bitrank->Save(file);
 337
 338     // Silly
 339     ulong *chars = new ulong[256/W+1];
 340     // Set up flags:
 341     for (ulong i = 0; i < 256; ++i)
 342         Tools::SetField(chars, 1, i, charMap->IsBitSet(i));
 343
 344     for (ulong offset = 0; offset < 256/W+1; ++offset)
 345     {
 346         if (std::fwrite(chars + offset, sizeof(ulong), 1, file) != 1)
 347             throw std::runtime_error("RLWaveletTree::Save(): file write error (chars).");
 348     }
 349
 350     delete [] chars;
 351 }
 352
 353 /**
 354  * Load from file
 355  */
 356 RLWaveletTree::RLWaveletTree(FILE *file)
 357 {
 358     if (std::fread(&(this->n), sizeof(ulong), 1, file) != 1)
 359         throw std::runtime_error("RLWaveletTree::Load(): file read error (n).");
 360     if (std::fread(&(this->maxlevel), sizeof(ulong), 1, file) != 1)
 361         throw std::runtime_error("RLWaveletTree::Load(): file read error (maxlevel).");
 362
 363     for (ulong offset = 0; offset < ((ulong)1 << maxlevel)+1; ++offset)
 364     {
 365         if (std::fread(this->leaves + offset, sizeof(ulong), 1, file) != 1)
 366             throw std::runtime_error("RLWaveletTree::Load(): file read error (leaves).");
 367     }
 368
 369     bitrank = new GapEncode(file);
 370
 371     // Silly
 372     ulong *chars = new ulong[256/W+1];
 373     // Read flags
 374     for (ulong offset = 0; offset < 256/W+1; ++offset)
 375     {
 376         if (std::fread(chars + offset, sizeof(ulong), 1, file) != 1)
 377             throw std::runtime_error("RLWaveletTree::Load(): file read error (chars).");
 378     }
 379     charMap = new BitRank(chars, 256, true);
 380
 381     prevRunLen = 0; // Init "buffer" values
 382     prevChar = 0;
 383     prevPos = 0;
 384 }