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