1 /***************************************************************************
2 * Copyright (C) 2007 by Veli Mäkinen *
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. *
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. *
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 ***************************************************************************/
21 #include "RLWaveletTree.h"
25 void RLWaveletTree::remapAlphabet(uchar *text, ulong n)
27 ulong *chars = new ulong[256/W+1];
29 for (i = 0; i < 256/W+1; ++i)
32 for (i = 0; i < n; ++i)
33 Tools::SetField(chars, 1, text[i], 1);
34 charMap = new BitRank(chars, 256, true);
37 for (i = 0; i < n; ++i)
38 text[i] = charMap->rank(text[i])-1;
41 RLWaveletTree::RLWaveletTree(uchar *D, ulong n) {
43 // std::cout << "Number of chars: " << charMap->rank(255) << ", maxlevel = " << Tools::CeilLog2(charMap->rank(255)) << std::endl;
46 this->maxlevel = Tools::CeilLog2(charMap->rank(255));
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 ++)
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++)
61 // compute the sizes of buckets
64 for (k=0; k<=maxlevel; k++)
65 bucket[k][value >> (maxlevel-k)]++;
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];
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];
83 // compute wtree using the information about boundaries
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);
90 // deallocate boundaries
91 for (k=0; k<=maxlevel; k++)
95 bitrank = new GapEncode(wtree, n*maxlevel, true);
101 RLWaveletTree::~RLWaveletTree() {
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;
112 ulong left = 0, right = (ulong)1 << maxlevel;
114 c = charMap->rank(c)-1;
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;
121 else { // go to the left
122 j -= bitrank->rank(j);
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;
132 else { // go to the left
133 j -= bitrank->rank(n*k+leaves[left]+j)-before;
134 right = (right+left)>>1;
136 // can never happen: if (left > right) return 0;
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
148 before = bitrank->rank(n*k+leaves[left]-1);
149 i = bitrank->select(before+i+1)-n*k-leaves[left];
151 else { // left side of parent
153 before = n*k+leaves[left]-bitrank->rank(n*k+leaves[left]-1);
154 i = bitrank->select0(before+i+1)-n*k-leaves[left];
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);
164 bool RLWaveletTree::IsCharAtPos(ulong c, ulong i) {
165 return charAtPos(i)==c;
169 ulong RLWaveletTree::charAtPos(ulong i)
171 // if (i >= n) // FIXME Replace with assert?
172 // std::cout << "RLWT::charAtPos parameter i out of range" << std::endl;
175 ulong left = 0, right = (ulong)1 << maxlevel;
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;
185 else { // go to the left
187 j -= rank_tmp; // bitrank->rank(j);
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;
198 else { // go to the left
200 j -= rank_tmp-before; //bitrank->rank(n*k+leaves[left]+j)-before;
201 right = (right+left)>>1;
203 // can never happen: if (left > right) return 0;
205 return charMap->select(c+1);
208 /* returns also the rank-1 of char at given position i */
209 ulong RLWaveletTree::charAtPos(ulong i, ulong *rank)
211 // if (i >= n) // FIXME Replace with assert?
212 // std::cout << "RLWT::charAtPos2 parameter i out of range" << std::endl;
215 ulong left = 0, right = (ulong)1 << maxlevel;
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;
225 else { // go to the left
227 j -= rank_tmp; //bitrank->rank(j);
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;
238 else { // go to the left
240 j -= rank_tmp-before; //bitrank->rank(n*k+leaves[left]+j)-before;
241 right = (right+left)>>1;
243 // can never happen: if (left > right) return 0;
246 return charMap->select(c+1);
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)
253 // Check for buffered data from previous call
254 if (i > prevPos && i <= prevPos + prevRunLen)
256 *rank = prevRank + (i-prevPos);
261 // if (i >= n) // FIXME Replace with assert?
262 // std::cout << "RLWT::charAtPos2 parameter i out of range" << std::endl;
265 ulong left = 0, right = (ulong)1 << maxlevel;
270 // first round separately to avoid extra "if"
271 if (bitrank->IsBitSet(j, &rank_tmp, &prevRunLen)) {
272 if (j+prevRunLen > n)
274 c = c | ((ulong)1 << (maxlevel-1));
275 j = rank_tmp-1; //bitrank->rank(j)-1;
278 else { // go to the left
280 if (j+prevRunLen > n)
282 j -= rank_tmp; //bitrank->rank(j);
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;
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;
304 j -= rank_tmp-before; //bitrank->rank(n*k+leaves[left]+j)-before;
305 right = (right+left)>>1;
307 // can never happen: if (left > right) return 0;
311 prevChar = charMap->select(c+1);
312 return charMap->select(c+1);
316 * Saving data fields:
319 ulong *leaves; // stores the leaf positions of characters 0..2^maxlevel-1
323 void RLWaveletTree::Save(FILE *file) const
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).");
330 for (ulong offset = 0; offset < ((ulong)1 << maxlevel)+1; ++offset)
332 if (std::fwrite(this->leaves + offset, sizeof(ulong), 1, file) != 1)
333 throw std::runtime_error("RLWaveletTree::Save(): file write error (leaves).");
339 ulong *chars = new ulong[256/W+1];
341 for (ulong i = 0; i < 256; ++i)
342 Tools::SetField(chars, 1, i, charMap->IsBitSet(i));
344 for (ulong offset = 0; offset < 256/W+1; ++offset)
346 if (std::fwrite(chars + offset, sizeof(ulong), 1, file) != 1)
347 throw std::runtime_error("RLWaveletTree::Save(): file write error (chars).");
356 RLWaveletTree::RLWaveletTree(FILE *file)
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).");
363 for (ulong offset = 0; offset < ((ulong)1 << maxlevel)+1; ++offset)
365 if (std::fread(this->leaves + offset, sizeof(ulong), 1, file) != 1)
366 throw std::runtime_error("RLWaveletTree::Load(): file read error (leaves).");
369 bitrank = new GapEncode(file);
372 ulong *chars = new ulong[256/W+1];
374 for (ulong offset = 0; offset < 256/W+1; ++offset)
376 if (std::fread(chars + offset, sizeof(ulong), 1, file) != 1)
377 throw std::runtime_error("RLWaveletTree::Load(): file read error (chars).");
379 charMap = new BitRank(chars, 256, true);
381 prevRunLen = 0; // Init "buffer" values