GapEncode.cpp

   1 /**
   2  * Gap encoding scheme
   3  * Niko Välimäki
   4  *
   5  */
   6
   7 #include "GapEncode.h"
   8 #include <stdexcept>
   9 using namespace std;
  10
  11
  12 GapEncode::GapEncode(ulong *B, ulong u, bool freeB = false)
  13 {
  14     ulong i, lastIndex, *bitRun, *runLength;
  15     this->u = u;
  16     this->numberOfRuns = 0;
  17     bitRun = new ulong[u/W + 1];
  18     for (i = 0; i < u/W + 1; i ++)
  19         bitRun[i]  = 0;
  20
  21     // Collect the runs of 1- or 0-bits
  22     totalLength = 0;
  23     uchar curBit = 0xFF;
  24     lastIndex = 0;
  25     for (i = 0; i < u; i ++)
  26         if (Tools::GetField(B, 1, i) != curBit)
  27         {
  28             if (curBit == 1)
  29                 totalLength += i - lastIndex;
  30             else if (Tools::GetField(B, 1, i))
  31                 Tools::SetField(bitRun, 1, i, 1);
  32
  33             lastIndex = i;
  34             curBit = Tools::GetField(B, 1, i);
  35         }
  36     // Handle the last run
  37     if (curBit == 1)
  38         totalLength += i - lastIndex;
  39
  40     runLength = new ulong[totalLength/W + 1];
  41     for (i = 0; i < totalLength/W + 1; i ++)
  42         runLength[i] = 0;
  43
  44     // Calculate run lengths
  45     totalLength = 0;
  46     curBit = 0xFF;
  47     lastIndex = 0;
  48     for (i = 0; i < u; i ++)
  49         if (Tools::GetField(B, 1, i) != curBit)
  50         {
  51             if (curBit == 1)
  52             {
  53                 totalLength += i - lastIndex;
  54                 Tools::SetField(runLength, 1, totalLength - 1, 1);
  55             }
  56             lastIndex = i;
  57             curBit = Tools::GetField(B, 1, i);
  58         }
  59     // Handle the last run
  60     if (curBit == 1)
  61     {
  62         totalLength += i - lastIndex;
  63         Tools::SetField(runLength, 1, totalLength - 1, 1);
  64     }
  65
  66     // Construct BSGAP structures
  67     this->bitRun = new BSGAP(bitRun, u, true, 0);
  68     numberOfRuns = this->bitRun->rank(u-1);
  69     this->runLength = new BSGAP(runLength, totalLength, true, 0);
  70
  71     if (freeB)
  72         delete [] B;
  73 }
  74
  75 GapEncode::~GapEncode()
  76 {
  77     // Free BSGAP structures
  78     delete bitRun;
  79     delete runLength;
  80 }
  81
  82 ulong GapEncode::rank(ulong i)
  83 {
  84     ulong pred, r = bitRun->rank(i);
  85     if (r == 0)
  86         return 0;
  87     pred = bitRun->select(r);
  88
  89     ulong bitsBeforeRun = 0;
  90     if (r > 1)
  91         bitsBeforeRun = runLength->select(r - 1) + 1;
  92
  93     ulong bitsInRun = runLength->select(r) + 1 - bitsBeforeRun;
  94     if (i - pred >= bitsInRun)
  95         return bitsBeforeRun + bitsInRun;
  96
  97     return bitsBeforeRun + i - pred + 1;
  98 }
  99
 100 ulong GapEncode::rank0(ulong i)
 101 {
 102     return i - rank(i) + 1;
 103 }
 104
 105 ulong GapEncode::select(ulong i)
 106 {
 107     throw logic_error("GapEncode::select() is not implemented");
 108     return 0;
 109 }
 110
 111 ulong GapEncode::select0(ulong i)
 112 {
 113     throw logic_error("GapEncode::select0() is not implemented");
 114     return 0; // To be implemented...
 115 }
 116
 117 /**
 118  * Saving data fields:
 119     ulong u;                // Universe size, number of bits in B
 120     ulong numberOfRuns;     // Number of 1-bit runs
 121     BSGAP *bitRun,          // Contains 1-bit at the start of each 1-bit-run
 122           *runLength;       // Contains run lengths for each 1-bit-run
 123     ulong totalLength;      // Total length of runs
 124 */
 125 void GapEncode::Save(FILE *file) const
 126 {
 127     if (std::fwrite(&(this->u), sizeof(ulong), 1, file) != 1)
 128         throw std::runtime_error("GapEncode::Save(): file write error (u).");
 129     if (std::fwrite(&(this->numberOfRuns), sizeof(ulong), 1, file) != 1)
 130         throw std::runtime_error("GapEncode::Save(): file write error (numberOfRuns).");
 131
 132     bitRun->Save(file);
 133     runLength->Save(file);
 134
 135     if (std::fwrite(&(this->totalLength), sizeof(ulong), 1, file) != 1)
 136         throw std::runtime_error("GapEncode::Save(): file write error (totalLength).");
 137 }
 138
 139 GapEncode::GapEncode(FILE *file)
 140 {
 141     if (std::fread(&(this->u), sizeof(ulong), 1, file) != 1)
 142         throw std::runtime_error("GapEncode::Load(): file read error (u).");
 143     if (std::fread(&(this->numberOfRuns), sizeof(ulong), 1, file) != 1)
 144         throw std::runtime_error("GapEncode::Load(): file read error (numberOfRuns).");
 145
 146     bitRun = new BSGAP(file);
 147     runLength = new BSGAP(file);
 148
 149     if (std::fread(&(this->totalLength), sizeof(ulong), 1, file) != 1)
 150         throw std::runtime_error("GapEncode::Load(): file read error (totalLength).");
 151 }