BSGAP.cpp

   1 /**
   2  * Binary-searchable gap encoding scheme (BSGAP)
   3  * Niko Välimäki
   4  *
   5  * Compile: g++ -Wall -o testBSGAP BSGAP.cpp Tools.cpp
   6  *
   7  * Based on:
   8  *
   9  * Ankur Gupta and Wing-Kai Hon and Rahul Shah and Jeffrey Scott Vitter. Compressed data structures:
  10  * Dictionaries and data-aware measures, Theor. Comput. Sci., Volume 387, Issue 3 (November 2007).
  11  */
  12
  13 #include "BSGAP.h"
  14 #include <stdexcept>
  15 #include <cassert>
  16 using namespace std;
  17
  18
  19 const uchar BSGAP::bit_table[] = {
  20     0,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,0,
  21     1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,6,0,1,0,2,0,1,0,3,0,1,0,
  22     2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,
  23     1,0,2,0,1,0,3,0,1,0,2,0,1,0,7,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,
  24     3,0,1,0,2,0,1,0,5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,
  25     1,0,6,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,
  26     2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1};
  27
  28 /**
  29  * samplerate == number of keys in each gap encoded binary search tree
  30  * default value is log^2 n where n is the number of keys.
  31  */
  32 BSGAP::BSGAP(ulong *B, ulong u_, bool freeB, ulong sampleRate)
  33     : u(u_), n(0), topCount(0), bitsInP(0), P(0), b(0), offsetA(0), firstKeyA(0)
  34 {
  35     // Count keys
  36     for (ulong i = 0; i < u; i ++)
  37         if (Tools::GetField(B, 1, i))
  38             this->n ++;
  39
  40     if (n == 0) // Sanity check
  41     {
  42         if (freeB)
  43             delete [] B;
  44
  45         return;
  46     }
  47
  48     unsigned log2n = Tools::FloorLog2(this->n);
  49     if (log2n < 2)
  50         log2n = 2;
  51     if (sampleRate == 0)
  52         this->b = log2n * log2n;
  53     else
  54         this->b = sampleRate;
  55
  56     ulong firstKey = 0, lastKey;
  57     // Find the first key
  58     while (!Tools::GetField(B, 1, firstKey))
  59         firstKey ++;
  60
  61     // Temp array of BSGAP structures
  62     ulong **tempP = new ulong * [n / b + 1];
  63     ulong *bitsInBSD = new ulong [n / b + 1];
  64     ulong totalBits = 0;
  65     // Number of nodes:
  66     this->topCount = n%b ? n/b + 1: n/b;
  67     this->firstKeyA = new BlockArray(topCount, Tools::CeilLog2(this->u));
  68
  69     for (ulong i = 0; i < n/b; i ++)
  70     {
  71         // Find the last key
  72         lastKey = firstKey;
  73         for (ulong k = 0; k < b + 1 && lastKey < u; lastKey ++)
  74             if (Tools::GetField(B, 1, lastKey))
  75                 k ++;
  76         lastKey --;
  77         // Sanity check
  78         if (lastKey > u)
  79         {
  80             cout << "error: lastKey = " << lastKey << ", u = " << u << endl;
  81             exit(0);
  82         }
  83
  84         // lastKey of the last substructure is this->u
  85         if (i == n/b - 1 && n%b == 0)
  86             lastKey = u;
  87
  88         tempP[i] = GetSubtree(B, firstKey, firstKey, lastKey, b, true, bitsInBSD[i]);
  89         totalBits += bitsInBSD[i];
  90         (*firstKeyA)[i] = firstKey;
  91
  92         firstKey = lastKey;
  93     }
  94
  95     if (n - (n/b) * b != 0)
  96     {
  97         // Find the first key
  98         while (!Tools::GetField(B, 1, firstKey))
  99             firstKey ++;
 100
 101         tempP[n/b] = GetSubtree(B, firstKey, firstKey, u, n - (n/b) * b, true, bitsInBSD[n/b]);
 102         totalBits += bitsInBSD[n/b];
 103         (*firstKeyA)[n/b] = firstKey;
 104     }
 105
 106     // Catenate binary trees into one bitvector
 107     this->P = new ulong[totalBits/W + 1];
 108     for (ulong i = 0; i < totalBits/W + 1; ++i)
 109         P[i] = 0;
 110     this->bitsInP = totalBits;
 111     ulong offset = 0;
 112     for (ulong i = 0; i < this->topCount ; ++i)
 113     {
 114         BitCopy(P, tempP[i], offset, bitsInBSD[i]);
 115         delete [] tempP[i];
 116         offset += bitsInBSD[i];
 117     }
 118     delete [] tempP;
 119
 120     // Pointers to bit vector P  (binary tree starting offsets)
 121     // FIXME Use more succinct representation?
 122     offset = 0;
 123     this->offsetA = new BlockArray(this->topCount, Tools::CeilLog2(totalBits));
 124     for (ulong i = 0; i < this->topCount ; ++i)
 125     {
 126         (*offsetA)[i] = offset;
 127         offset += bitsInBSD[i];
 128     }
 129     delete [] bitsInBSD;
 130
 131     if (freeB)
 132         delete [] B;
 133 }
 134
 135 BSGAP::~BSGAP()
 136 {
 137         delete [] P;
 138         delete offsetA;
 139         delete firstKeyA;
 140 }
 141
 142 ulong BSGAP::rank(ulong i)
 143 {
 144     if (n == 0) // Trivial case
 145         return 0;
 146     if (i >= u)
 147         throw std::out_of_range("BSGAP::rank(): Parameter i was out of range");
 148     if ((*firstKeyA)[0] > i)
 149         return 0;
 150     ulong l = 0, r = topCount - 1;
 151     while (l < r)   // Binary search on the array firstKeyA
 152     {
 153         ulong mid = l + (r - l)/2lu;
 154         if ((*firstKeyA)[mid] < i)
 155             l = mid + 1;
 156         else
 157             r = mid;
 158     }
 159     if ((*firstKeyA)[l] > i && l == 0)
 160         return 0;
 161     if ((*firstKeyA)[l] == i)
 162         return 1 + l * b;
 163     if ((*firstKeyA)[l] > i)
 164         l --;
 165
 166     ulong lastKey, firstKey = (*firstKeyA)[l];
 167     if (l == topCount - 1)
 168         lastKey = u;
 169     else
 170         lastKey = (*firstKeyA)[l + 1];
 171     ulong nSub = b;
 172     if (l == topCount - 1 && n%b)
 173         nSub = n%b;
 174
 175     ulong k = rank(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true);
 176     // Add rank of substructures P[0], P[1], ..., P[j - 2]
 177     return l * b + k;
 178 }
 179
 180 ulong BSGAP::rank0(ulong i)
 181 {
 182     return i - rank(i) + 1;
 183 }
 184
 185 bool BSGAP::IsBitSet(ulong i)
 186 {
 187     if (n == 0) // Trivial case
 188         return false;
 189     if (i >= u)
 190         return false; // FIXME throw std::out_of_range("BSGAP::rank(): Parameter i was out of range");
 191     if ((*firstKeyA)[0] > i)
 192         return false;
 193     ulong l = 0, r = topCount - 1;
 194     while (l < r)   // Binary search on the array firstKeyA
 195     {
 196         ulong mid = l + (r - l)/2lu;
 197         if ((*firstKeyA)[mid] < i)
 198             l = mid + 1;
 199         else
 200             r = mid;
 201     }
 202     if ((*firstKeyA)[l] > i && l == 0)
 203         return false;
 204     if ((*firstKeyA)[l] == i)
 205         return true;
 206     if ((*firstKeyA)[l] > i)
 207         l --;
 208
 209     ulong lastKey, firstKey = (*firstKeyA)[l];
 210     if (l == topCount - 1)
 211         lastKey = u;
 212     else
 213         lastKey = (*firstKeyA)[l + 1];
 214     ulong nSub = b;
 215     if (l == topCount - 1 && n%b)
 216         nSub = n%b;
 217
 218     return IsBitSet(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true);
 219 }
 220
 221 /**
 222  * Returns also rank_1(i) via the second parameter.
 223  */
 224 bool BSGAP::IsBitSet(ulong i, ulong *resultRank)
 225 {
 226     *resultRank = 0;
 227     if (n == 0) // Trivial case
 228         return false;
 229     if (i >= u)
 230     {
 231         *resultRank = n;
 232         return false; // FIXME throw std::out_of_range("BSGAP::rank(): Parameter i was out of range");
 233     }
 234     if ((*firstKeyA)[0] > i)
 235         return false;
 236     ulong l = 0, r = topCount - 1;
 237     while (l < r)   // Binary search on the array firstKeyA
 238     {
 239         ulong mid = l + (r - l)/2lu;
 240         if ((*firstKeyA)[mid] < i)
 241             l = mid + 1;
 242         else
 243             r = mid;
 244     }
 245     if ((*firstKeyA)[l] > i && l == 0)
 246         return false;
 247     if ((*firstKeyA)[l] == i)
 248     {
 249         *resultRank =  1 + l * b;
 250         return true;
 251     }
 252     if ((*firstKeyA)[l] > i)
 253         l --;
 254
 255     ulong lastKey, firstKey = (*firstKeyA)[l];
 256     if (l == topCount - 1)
 257         lastKey = u;
 258     else
 259         lastKey = (*firstKeyA)[l + 1];
 260     ulong nSub = b;
 261     if (l == topCount - 1 && n%b)
 262         nSub = n%b;
 263
 264     bool result = IsBitSet(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true, resultRank);
 265     *resultRank += l * b;
 266     return result;
 267 }
 268
 269
 270 ulong BSGAP::select(ulong i)
 271 {
 272     if (n == 0)
 273         return 0;
 274     if (i == 0)
 275         return 0;
 276     if (i > n)
 277         throw std::out_of_range("BSGAP::select(): Parameter i was out of range");
 278
 279     ulong j = i/b;
 280     if (i%b == 0)
 281         j --;
 282     ulong lastKey, firstKey = (*firstKeyA)[j];
 283     if (j == topCount - 1)
 284         lastKey = u;
 285     else
 286         lastKey = (*firstKeyA)[j + 1];
 287     ulong nSub = b;
 288     if (j == topCount - 1 && n%b)
 289         nSub = n%b;
 290
 291     return select(P, (*offsetA)[j], firstKey, i - j * b, firstKey, lastKey, nSub, true);
 292 }
 293
 294 ulong BSGAP::select0(ulong i)
 295 {
 296     if (n == 0) // Trivial case
 297         return i-1;
 298     if (i > u - n)
 299         throw std::out_of_range("BSGAP::select0(): Parameter i was out of range");
 300
 301     ulong l = 0, r = topCount - 1;
 302     while (l < r)   // Binary search on the array firstKeyA
 303     {
 304         ulong mid = l + (r - l)/2lu;
 305         if ((*firstKeyA)[mid] - mid*b < i) // Number of zeros before [mid]
 306             l = mid + 1;
 307         else
 308             r = mid;
 309     }
 310     if (l == 0 && (*firstKeyA)[l] >= i)
 311         return i-1;
 312     if (l > 0 && (*firstKeyA)[l] - l*b >= i)
 313         --l; // FIXME: TEST
 314
 315     i -= (*firstKeyA)[l] - l*b;
 316
 317     ulong lastKey, firstKey = (*firstKeyA)[l];
 318     if (l == topCount - 1)
 319         lastKey = u;
 320     else
 321         lastKey = (*firstKeyA)[l + 1];
 322     ulong nSub = b;
 323     if (l == topCount - 1 && n%b)
 324         nSub = n%b;
 325
 326     return select0(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true);
 327 }
 328
 329 ulong BSGAP::select0(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 330 {
 331     while (1)
 332     {
 333         if (n == 0)
 334             return l+i;
 335
 336         // Check if subtree is full  --- FIXME Make this test an assert()!
 337         if (IsSubtreeFull(firstKey, l, r, n))
 338             throw runtime_error("BSGAP::select0(): subtree was full"); // Should not happen // (l == firstKey ? l + i - 1: l + i);
 339
 340         // Parse key value
 341         ulong key;
 342         ulong x = DeltaDecode(B, offset);
 343         if (leftChild)
 344             key = r - x;
 345         else
 346             key = l + x;
 347
 348         if (numberOfZeros(firstKey, l, key, n) < i)
 349         {
 350             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 351             i = i - numberOfZeros(firstKey, l, key, n);
 352             l = key;
 353             n = n - n/2 - 1;
 354             leftChild = false;
 355         }
 356         else
 357         {
 358             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 359             r = key;
 360             n = n/2;
 361             leftChild = true;
 362         }
 363     }
 364 }
 365
 366
 367 ulong BSGAP::rank(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 368 {
 369     // Number of keys in subtrees on left-side
 370     ulong result = 0;
 371
 372     while (1)
 373     {
 374
 375         if (n == 0)
 376             return result;
 377
 378         // Check if subtree is full
 379         if (IsSubtreeFull(firstKey, l, r, n))
 380             return result + (l == firstKey ? i - l + 1 : i - l);
 381
 382         // Parse key value
 383         ulong key;
 384         ulong x = DeltaDecode(B, offset);
 385         if (leftChild)
 386             key = r - x;
 387         else
 388             key = l + x;
 389         if (key == i)
 390             return result + n/2 + 1;
 391
 392         if (key < i)
 393         {
 394             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 395             result += n/2 + 1;
 396             l = key;
 397             n = n - n/2 - 1;
 398             leftChild = false;
 399         }
 400         else
 401         {
 402             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 403             r = key;
 404             n = n/2;
 405             leftChild = true;
 406         }
 407     }
 408 }
 409
 410
 411 bool BSGAP::IsBitSet(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 412 {
 413     // Number of keys in subtrees on left-side
 414     ulong result = 0;
 415
 416     while (1)
 417     {
 418
 419         if (n == 0)
 420             return false;
 421
 422         // Check if subtree is full
 423         if (IsSubtreeFull(firstKey, l, r, n))
 424             return result + (l == firstKey ? i - l + 1 : i - l);
 425
 426         // Parse key value
 427         ulong key;
 428         ulong x = DeltaDecode(B, offset);
 429         if (leftChild)
 430             key = r - x;
 431         else
 432             key = l + x;
 433         if (key == i)
 434             return true;
 435
 436         if (key < i)
 437         {
 438             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 439             result += n/2 + 1;
 440             l = key;
 441             n = n - n/2 - 1;
 442             leftChild = false;
 443         }
 444         else
 445         {
 446             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 447             r = key;
 448             n = n/2;
 449             leftChild = true;
 450         }
 451     }
 452 }
 453
 454 // Returns also the rank of i
 455 bool BSGAP::IsBitSet(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild, ulong *resultRank)
 456 {
 457     // Number of keys in subtrees on left-side
 458     ulong result = 0;
 459
 460     while (1)
 461     {
 462
 463         if (n == 0)
 464         {
 465             *resultRank = result;
 466             return false;
 467         }
 468
 469         // Check if subtree is full
 470         if (IsSubtreeFull(firstKey, l, r, n))
 471         {
 472             *resultRank = result + (l == firstKey ? i - l + 1 : i - l);
 473             return true;
 474         }
 475
 476         // Parse key value
 477         ulong key;
 478         ulong x = DeltaDecode(B, offset);
 479         if (leftChild)
 480             key = r - x;
 481         else
 482             key = l + x;
 483         if (key == i)
 484         {
 485             *resultRank = result + n/2 + 1;
 486             return true;
 487         }
 488
 489         if (key < i)
 490         {
 491             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 492             result += n/2 + 1;
 493             l = key;
 494             n = n - n/2 - 1;
 495             leftChild = false;
 496         }
 497         else
 498         {
 499             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 500             r = key;
 501             n = n/2;
 502             leftChild = true;
 503         }
 504     }
 505 }
 506
 507
 508 ulong BSGAP::select(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 509 {
 510     while (1)
 511     {
 512         if (n == 0)
 513             return 0;
 514
 515         // Check if subtree is full
 516         if (IsSubtreeFull(firstKey, l, r, n))
 517             return (l == firstKey ? l + i - 1: l + i);
 518
 519         // Parse key value
 520         ulong key;
 521         ulong x = DeltaDecode(B, offset);
 522         if (leftChild)
 523             key = r - x;
 524         else
 525             key = l + x;
 526         if (n/2 + 1 == i)
 527             return key;
 528         if (n/2 + 1 < i)
 529         {
 530             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 531             i = i - n/2 - 1;
 532             l = key;
 533             n = n - n/2 - 1;
 534             leftChild = false;
 535         }
 536         else
 537         {
 538             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 539             r =  key;
 540             n = n/2;
 541             leftChild = true;
 542         }
 543     }
 544 }
 545
 546 ulong * BSGAP::GetSubtree(ulong *B, ulong firstKey, ulong l, ulong r, ulong n, bool leftChild, ulong &bits)
 547 {
 548     bits = 0;
 549     if (n == 0)
 550         return 0;
 551
 552     // Check if subtree is full
 553     if (IsSubtreeFull(firstKey, l, r, n))
 554         return 0;
 555
 556     // Pick the median
 557     ulong key = FindMedian(B, firstKey, l, r, n);
 558
 559     ulong *leftSub, *rightSub, leftBits, rightBits;
 560     leftSub = GetSubtree(B, firstKey, l, key, n/2, true, leftBits);
 561     rightSub = GetSubtree(B, firstKey, key, r, n - n/2 - 1, false, rightBits);
 562
 563     // Encode key value
 564     ulong *keyDelta, keyBits;
 565     if (leftChild)
 566         keyDelta = DeltaEncode(r - key, keyBits, true);
 567     else
 568         keyDelta = DeltaEncode(key - l, keyBits, true);
 569
 570     // Encode jump offset if left and right subtrees exists
 571     ulong *jumpOffset = 0, jumpBits = 0;
 572     if (leftBits != 0 && rightBits != 0)
 573         jumpOffset = DeltaEncode(leftBits, jumpBits);
 574
 575     // bits is the sum of keyBits, jumpBits, leftBits and rightBits
 576     bits = keyBits + jumpBits + leftBits + rightBits;
 577     ulong *output = new ulong[bits / W + 1];
 578     for (ulong i = 0; i < bits/W+1; ++i)
 579         output[i] = 0;
 580
 581     /**
 582      * Write "output"
 583      */
 584     BitCopy(output, keyDelta, 0, keyBits);
 585     delete [] keyDelta;
 586     ulong offset = keyBits;
 587
 588     if (jumpBits != 0)
 589     {
 590         BitCopy(output, jumpOffset, offset, jumpBits);
 591         offset += jumpBits;
 592         delete [] jumpOffset;
 593     }
 594
 595     BitCopy(output, leftSub, offset, leftBits);
 596     offset += leftBits;
 597     BitCopy(output, rightSub, offset, rightBits);
 598     offset += rightBits;
 599
 600     assert(offset == bits);
 601     if (leftBits != 0)
 602         delete [] leftSub;
 603     if (rightBits != 0)
 604         delete [] rightSub;
 605     return output;
 606 }
 607
 608 void BSGAP::BitCopy(ulong *dest, ulong *src, ulong offset, ulong len)
 609 {
 610     if (len == 0)
 611         return;
 612     ulong i;
 613
 614     if (len <= W)
 615     {
 616         Tools::SetVariableField(dest, len, offset, *src);
 617         return;
 618     }
 619
 620     for (i = 0; i < len/W; i ++)
 621         Tools::SetVariableField(dest, W, offset + i * W, src[i]);
 622
 623     if (i * W < len)
 624         Tools::SetVariableField(dest, len - i*W, offset + i * W, src[i]);
 625 }
 626
 627 ulong BSGAP::FindMedian(ulong *B, ulong firstKey, ulong l, ulong r, ulong n)
 628 {
 629     // Linear scan: slow but affects only the construction time
 630     n = n/2 + 1;
 631     if (l != firstKey)
 632         l ++;   // Skip left boundary
 633
 634     while (n && l <= r)
 635     {
 636         if (Tools::GetField(B, 1, l))
 637             n --;
 638         l ++;
 639     }
 640     return l - 1;
 641 }
 642
 643 ulong BSGAP::DeltaDecode(ulong *B, ulong &offset)
 644 {
 645     // Dense coding
 646     // http://www.dcc.uchile.cl/~gnavarro/ps/ir06.pdf
 647     const unsigned s = 180; // FIXME These should be class variables
 648     const unsigned c = 256 - s; // FIXME Choose <s,c> values!
 649     const unsigned b = 8; // one byte.
 650
 651     ulong i = 0;
 652     ulong base = 0;
 653     ulong tot = s;
 654     while (Tools::GetVariableField(B, b, offset) < c)
 655     {
 656         i = i*c + Tools::GetVariableField(B, b, offset);
 657         offset += b;
 658         base += tot;
 659         tot = tot * c;
 660     }
 661     i = i * s + (Tools::GetVariableField(B, b, offset) - c);
 662     offset += b;
 663     return i + base;
 664 }
 665
 666
 667 ulong * BSGAP::DeltaEncode(ulong value, ulong &bits, bool onlyPositive, bool negative)
 668 {
 669     // Dense coding
 670     // http://www.dcc.uchile.cl/~gnavarro/ps/ir06.pdf
 671     const unsigned s = 180;  // FIXME These should be class variables
 672     const unsigned c = 256 - s; // FIXME Choose <s,c> values!
 673     const unsigned b = 8; // one byte
 674
 675     // Calculate size first:
 676     bits = b;
 677     ulong x = value / s;
 678     while (x > 0)
 679     {
 680         --x;
 681         bits += b;
 682         x = x/c;
 683     }
 684
 685     // bits is now the length of the whole codeword
 686     ulong *B = new ulong[bits/W + 1];
 687
 688     // output codeword (backwards):
 689     unsigned offset = bits - b;
 690     ulong output = c + (value % s);
 691     Tools::SetVariableField(B, b, offset, output);
 692
 693     x = value / s;
 694     while (x > 0)
 695     {
 696         --x;
 697         offset -= b;
 698
 699         output = x % c;
 700         Tools::SetVariableField(B, b, offset, output);
 701
 702         x = x/c;
 703     }
 704
 705     return B;
 706 }
 707
 708
 709 /**
 710  * Saving the following data fields:
 711  *   ulong u, n;             // Universe size, number of 1-bits in B
 712  *   ulong topCount;   // Top structure and the number of substructures
 713  *   ulong bitsInP;
 714  *   ulong *P;               // Pointer to BSGAP structures
 715  *   unsigned b;             // Subdictionary size (\log^2 n)
 716  *   BlockArray *offsetA;    // Array of pointers (into bitvector P)
 717  *   BlockArray *firstKeyA;  // Array of first key positions of the substructures
 718  */
 719 void BSGAP::Save(FILE *file) const
 720 {
 721     if (std::fwrite(&(this->u), sizeof(ulong), 1, file) != 1)
 722         throw std::runtime_error("BSGAP::Save(): file write error (u).");
 723
 724     if (std::fwrite(&(this->n), sizeof(ulong), 1, file) != 1)
 725         throw std::runtime_error("BSGAP::Save(): file write error (n).");
 726
 727     if (n == 0)
 728         return; // Done.
 729
 730     if (std::fwrite(&(this->topCount), sizeof(ulong), 1, file) != 1)
 731         throw std::runtime_error("BSGAP::Save(): file write error (topCount).");
 732
 733     if (std::fwrite(&(this->bitsInP), sizeof(ulong), 1, file) != 1)
 734         throw std::runtime_error("BSGAP::Save(): file write error (bitsInP).");
 735
 736     for (ulong offset = 0; offset < bitsInP/W+1; offset ++)
 737     {
 738         if (std::fwrite((this->P)+offset, sizeof(ulong), 1, file) != 1)
 739             throw std::runtime_error("BSGAP::Save(): file write error (P).");
 740             }
 741
 742     if (std::fwrite(&(this->b), sizeof(unsigned), 1, file) != 1)
 743         throw std::runtime_error("BSGAP::Save(): file write error (b).");
 744
 745     offsetA->Save(file);
 746     firstKeyA->Save(file);
 747 }
 748
 749 /**
 750  * Load from file
 751  */
 752 BSGAP::BSGAP(FILE *file)
 753     : u(0), n(0), topCount(0), bitsInP(0), P(0), b(0), offsetA(0), firstKeyA(0)
 754 {
 755     if (std::fread(&(this->u), sizeof(ulong), 1, file) != 1)
 756         throw std::runtime_error("BSGAP::Load(): file read error (u).");
 757
 758     if (std::fread(&(this->n), sizeof(ulong), 1, file) != 1)
 759         throw std::runtime_error("BSGAP::Load(): file read error (n).");
 760
 761     if (n == 0)
 762         return; // Done.
 763
 764     if (std::fread(&(this->topCount), sizeof(ulong), 1, file) != 1)
 765         throw std::runtime_error("BSGAP::Load(): file read error (topCount).");
 766
 767     if (std::fread(&(this->bitsInP), sizeof(ulong), 1, file) != 1)
 768         throw std::runtime_error("BSGAP::Load(): file read error (bitsInP).");
 769
 770     P = new ulong[bitsInP/W+1];
 771     for (ulong offset = 0; offset < bitsInP/W+1; offset ++)
 772     {
 773         if (std::fread(this->P+offset, sizeof(ulong), 1, file) != 1)
 774             throw std::runtime_error("BSGAP::Load(): file read error (P).");
 775     }
 776
 777     if (std::fread(&(this->b), sizeof(unsigned), 1, file) != 1)
 778         throw std::runtime_error("BSGAP::Load(): file read error (b).");
 779
 780     offsetA = new BlockArray(file);
 781     firstKeyA = new BlockArray(file);
 782 }
 783