LPAQ

Posted by Borelset on June 7, 2022

PAQ系列压缩器把数据压缩问题完全转换为了数据建模问题,根据数据建模进行逐个bit的预测,并结合算术编码来对数据进行压缩。 LPAQ是一个基于PAQ8l的简化版本。

Predictor

LPAQ的核心是Predictor。它一方面负责根据输入的数据来进行学习,调整数据建模;另一方面则根据过去的输入结果来对下一个bit进行预测。 Predictor类的成员如下

class Predictor {
    int pr;  // next prediction

    U8 t0[0x10000];  // order 1 cxt -> mode  // 保存order1的context对应的历史状态
    HashTable<16> t;  // cxt -> mode  //保存order 2, 3, 4, 6, unigram的context对应的历史状态。
    int c0 = 1;  // last 0-7 bits with leading 1 // 最新的若干bit
    int c4 = 0;  // last 4 bytes  // 最近的4字节
    U8 *cp[6] = {t0, t0, t0, t0, t0, t0};  // pointer to bit history // 分别对应order2,3,4,6,unigram的状态
    int bcount = 0;  // bit count
    StateMap sm[6];
    APM a1, a2;
    U32 h[6];
    Mixer m;
    MatchModel mm;  // predicts next bit by matching context

Predictor的原理其实就是“统计出现某种情况(数据片段)的时候下一个bit是什么,然后再次出现这种情况(数据片段)的时候,就根据过去的统计结果来进行预测”。 而需要针对性统计的数据片段在这里被称为context。

Predictor通过7种context来进行预测,分别是order 1, order 2, order 3, order 4, order 6, unigram和MatchModel。

  • order n这种context可以认为是一个长度为n字节的前缀。比如说,如果我有一个字符串“abcdefg”,那么当需要预测g的下一个bit时,order 1的context就是“g”,order 2的context就是“fg”,order 3的context就是“efg”,以此类推。
  • unigram种context是设计出来专门针对字母型语言的文本数据的,它会忽略非字母数据(ACSII),并且忽略文本的大小写,并且在遇到文本空格时截断。
  • MatchModel这种context会利用过去的数据流,如果遇到了过去出现过的重复片段,则会直接根据过去的历史来进行预测。

在上述的这些类成员中

  • t0被用来记录order 1 context的统计结果。
  • t被用来记录order 2, order 3, order 4, order 6 context的统计结果。
  • mm被用来记录MatchModel context的统计结果。
  • c0和c4记录了最近的若干bit和字节,为了方便计算得到context。(虽然我们说order n这种context就是最近的n个字节的前缀,但是在代码中会根据这个前缀计算出一个哈希值代替context本身)
  • cp是六个指针,分别对应order 1, order 2, order 3, order 4, order 6, unigram context所指向的记录。比如说,目前的order 1 context是“b”,有关于“b”的统计结果被保存于t0中的某个位置pos。那么此时cp[0]就应该是t0[pos]。
  • bcount用来记录当前处理的bit是一个字节的第几个bit。(Predictor逐bit的进行学习和预测)
  • sm用来进行实现从统计到预测概率的转换。一个数组中又饿6个StateMap,分别对应order 1, order 2, order 3, order 4, order 6和unigram。MatchModel的数据结构里有一个属于自己的StateMap,因此不在sm这个数组里。
  • a1和a2负责最终依次对预测概率的修正
  • h分别保存order 1, order 2, order 3, order 4, order 6, unigram context计算出的哈希值。
  • m负责把7种不同context得到的预测概率进行加权混合。

StateMap

StateMap负责根据统计结果来进行概率预测。 其代码如下所示

class StateMap {
protected:
    const int N;  // Number of contexts  // state的总数量
    int cxt;      // Context of last prediction // 当前的state
    U32 *t;       // cxt -> prediction in high 22 bits, count in low 10 bits // state对应到的概率和计数
    static int dt[1024];  // i -> 16K/(i+3)  // 学习率
    void update(int y, int limit) {
      assert(cxt >= 0 && cxt < N);
      int n = t[cxt] & 1023, p = t[cxt] >> 10;  // count, prediction // 前22bit代表概率,后10bit是计数。
      if (n < limit) ++t[cxt];   // 计数+1
      else t[cxt] = t[cxt] & 0xfffffc00 | limit;  // t[cxt]&0xfffffc00 表示提取前22bit,|limit表示把计数填充为limit值。
      t[cxt] += (((y << 22) - p) >> 3) * dt[n] &
                0xfffffc00; // 更新概率值 y<<22-p 表示预测误差,右移三位乘以dt[n]应该和学习率有关,然后&0xfffffc00表示只保留高22bit。
    }
public:
    StateMap(int n = 256);
    
    ~StateMap() {
      free(t);
    }
    
    // update bit y (0..1), predict next bit in context cx
    int p(int y, int cx, int limit = 1023) {
      assert(y >> 1 == 0);
      assert(cx >= 0 && cx < N);
      assert(limit > 0 && limit < 1024);
      update(y, limit);
      return t[cxt = cx] >> 20; // 右移20位,实际等于t[cxt]中的概率原本是22位。也右移了10位,输出12bit的概率。
    }
};

int StateMap::dt[1024] = {0};

StateMap::StateMap(int n) : N(n), cxt(0) {
  alloc(t, N);
  for (int i = 0; i < N; ++i)
    t[i] = 1 << 31;
  if (dt[0] == 0)
    for (int i = 0; i < 1024; ++i)
      dt[i] = 16384 / (i + i + 3); // 这个dt[i]应该是学习率?
}

StateMap类包含了四个成员。

  • N,它代表了这个StateMap所记录的state的数量
  • cxt,它代表了当前的state
  • t,它记录了某个state对应的概率,以及出现总次数
  • dt,它代表某个state在出现第k次时,修正对应概率时的学习率

我们先来看StateMap的构造函数

StateMap::StateMap(int n) : N(n), cxt(0) {
  alloc(t, N);
  for (int i = 0; i < N; ++i)
    t[i] = 1 << 31;
  if (dt[0] == 0)
    for (int i = 0; i < 1024; ++i)
      dt[i] = 16384 / (i + i + 3); // 这个dt[i]应该是学习率?
}

这个函数接受一个参数n。这个参数会直接赋值给N,代表了StateMap所记录的state的数量。

alloc是一个封装了malloc的函数,它会根据t的指针类型,来获取N个t对应类型的大小,并将地址保存在t。 比如,如果t是uint64_t的类型,那么这个函数就会获取sizeof(uint64_t)N的内存空间。 这里t是int*,所以最终会获取4N字节的内存空间。

给t获取内存空间后,于是就开始对t初始化。t是一个数组,这个数组中的每一个元素都对应了一个状态。 初始化的过程通过一个for循环,依次给数组t中的所有元素赋初值1«31(也就是0xc0000000)。 这里我们先来解释这个数值的含义。 我们知道数组t中的元素是int类型,也就是32bit。StateMap将这32bit划分为两个部分:

  • 前22bit:表示这个state对应的下一个bit是1的概率大小。
    • 最大值为二进制的22个1,也就是0x3fffff,对应下一个bit是1的概率为1.0;
    • 最小值为二进制的22个0,也就是0x000000,对应下一个bit是1的概率为0.0;
    • 中间值为二进制的21个1,也就是0x1fffff,对应下一个bit是1的概率为0.5;
  • 后10bit:表示这个state出现的次数 可以发现0xc0000000的前22个bit就是0x1fffff;而后10个bit就是0x000。 也就是说,0xc0000000这个数字的含义其实就是,当前状态对应的概率为0.5,并且这个状态以前从未出现过。

dt则表示对某个状态对应的概率进行更新的幅度有多大,这只是一个系数,没有具体的含义。 但可以发现,这个系数的设置是倾向于某个状态出现的次数越多,概率更新的幅度越小。

接下来我们再看StateMap生成对应概率的函数

int p(int y, int cx, int limit = 1023) {
  assert(y >> 1 == 0);
  assert(cx >= 0 && cx < N);
  assert(limit > 0 && limit < 1024);
  update(y, limit);
  return t[cxt = cx] >> 20; // 右移20位,实际等于t[cxt]中的概率原本是22位。也右移了10位,输出12bit的概率。
}

这个函数接受三个参数

  • 第一个参数y,用来接收当前处理到的bit
  • 第二个参数cx,用来接收当前的状态
  • 第三个参数limit,表示在数组t的元素中,记录状态出现次数的上限。我们前面介绍过,记录次数使用10个bit,所以这个上限就是1023。

这个函数首先会调用update

void update(int y, int limit) {
  assert(cxt >= 0 && cxt < N);
  int n = t[cxt] & 1023, p = t[cxt] >> 10;  // count, prediction // 前22bit代表概率,后10bit是计数。
  if (n < limit) ++t[cxt];   // 计数+1
  else t[cxt] = t[cxt] & 0xfffffc00 | limit;  // t[cxt]&0xfffffc00 表示提取前22bit,|limit表示把计数填充为limit值。
  // 更新概率值 y<<22-p 表示预测误差,右移三位乘以dt[n]应该和学习率有关,然后&0xfffffc00表示只保留高22bit。
  t[cxt] += (((y << 22) - p) >> 3) * dt[n] & 0xfffffc00; 
}

这个中的cxt应该还是上一次调用时留下来的,也就是上一个状态。而这里的y是当前正在处理的bit。 所以,这个函数其实就是在对上一个状态所作出的预测进行修正。

t[cxt]就是上一个状态所对应的元素。 n获取了t[cxt]的后10个bit,也就是该状态出现的次数; 而p则获取了t[cxt]的前22个bit,也就是该状态所对应的概率。

接下来就检查n是否超出了最大值,如果没有超出,那么++t[cxt]其实就是相当于在后10个bit上+1,记录该状态出现的次数+1; 如果超出了最大限制,那么就先获取t[cxt]前22bit的概率部分(t[cxt] & 0xfffffc00),然后再把一个出现次数的最大值放在后10个bit上(| limit)。实际上这个步骤感觉多余。

最后还需要更新t[cxt]的概率部分。 ((y«22)-p)其实就是上一个状态的预测误差,这里y左移22位是因为p是22bit的。 这里将预测误差右移3位(相当于除以8),乘以一个预设系数dt[n],再与(&)0xfffffc00,最终得到一个t[cxt]的概率部分的更新值。 0xfffffc00一方面保证了这个更新值不会影响到t[cxt]中记录出现次数的部分,另一方面也相当于让更新值再右移10位。 也就是说对于前面提到的p这个int型整数,它在这一行实际的更新值为(((y«22)-p) » 3) * dt[n] » 10。

接着我们回到int p(int y, int cx, int limit = 1023)这个函数。该函数最后一行会更新cxt,把它设置为新的一个状态,然后再获取新的状态下,下一个bit为1的预测概率。

HashTable

LPAQ自己实现了一个哈希表,来记录order 2, order 3, order 4, order 6,和unigram的context上次出现时的所对应的状态。 这里的状态指的是LPAQ预定义的253种状态(和PAQ8中的相同),如下方所示。

static const U8 State_table[256][2]={
{  1,  2},{  3,  5},{  4,  6},{  7, 10},{  8, 12},{  9, 13},{ 11, 14}, // 0
{ 15, 19},{ 16, 23},{ 17, 24},{ 18, 25},{ 20, 27},{ 21, 28},{ 22, 29}, // 7
{ 26, 30},{ 31, 33},{ 32, 35},{ 32, 35},{ 32, 35},{ 32, 35},{ 34, 37}, // 14
{ 34, 37},{ 34, 37},{ 34, 37},{ 34, 37},{ 34, 37},{ 36, 39},{ 36, 39}, // 21
{ 36, 39},{ 36, 39},{ 38, 40},{ 41, 43},{ 42, 45},{ 42, 45},{ 44, 47}, // 28
{ 44, 47},{ 46, 49},{ 46, 49},{ 48, 51},{ 48, 51},{ 50, 52},{ 53, 43}, // 35
{ 54, 57},{ 54, 57},{ 56, 59},{ 56, 59},{ 58, 61},{ 58, 61},{ 60, 63}, // 42
{ 60, 63},{ 62, 65},{ 62, 65},{ 50, 66},{ 67, 55},{ 68, 57},{ 68, 57}, // 49
{ 70, 73},{ 70, 73},{ 72, 75},{ 72, 75},{ 74, 77},{ 74, 77},{ 76, 79}, // 56
{ 76, 79},{ 62, 81},{ 62, 81},{ 64, 82},{ 83, 69},{ 84, 71},{ 84, 71}, // 63
{ 86, 73},{ 86, 73},{ 44, 59},{ 44, 59},{ 58, 61},{ 58, 61},{ 60, 49}, // 70
{ 60, 49},{ 76, 89},{ 76, 89},{ 78, 91},{ 78, 91},{ 80, 92},{ 93, 69}, // 77
{ 94, 87},{ 94, 87},{ 96, 45},{ 96, 45},{ 48, 99},{ 48, 99},{ 88,101}, // 84
{ 88,101},{ 80,102},{103, 69},{104, 87},{104, 87},{106, 57},{106, 57}, // 91
{ 62,109},{ 62,109},{ 88,111},{ 88,111},{ 80,112},{113, 85},{114, 87}, // 98
{114, 87},{116, 57},{116, 57},{ 62,119},{ 62,119},{ 88,121},{ 88,121}, // 105
{ 90,122},{123, 85},{124, 97},{124, 97},{126, 57},{126, 57},{ 62,129}, // 112
{ 62,129},{ 98,131},{ 98,131},{ 90,132},{133, 85},{134, 97},{134, 97}, // 119
{136, 57},{136, 57},{ 62,139},{ 62,139},{ 98,141},{ 98,141},{ 90,142}, // 126
{143, 95},{144, 97},{144, 97},{ 68, 57},{ 68, 57},{ 62, 81},{ 62, 81}, // 133
{ 98,147},{ 98,147},{100,148},{149, 95},{150,107},{150,107},{108,151}, // 140
{108,151},{100,152},{153, 95},{154,107},{108,155},{100,156},{157, 95}, // 147
{158,107},{108,159},{100,160},{161,105},{162,107},{108,163},{110,164}, // 154
{165,105},{166,117},{118,167},{110,168},{169,105},{170,117},{118,171}, // 161
{110,172},{173,105},{174,117},{118,175},{110,176},{177,105},{178,117}, // 168
{118,179},{110,180},{181,115},{182,117},{118,183},{120,184},{185,115}, // 175
{186,127},{128,187},{120,188},{189,115},{190,127},{128,191},{120,192}, // 182
{193,115},{194,127},{128,195},{120,196},{197,115},{198,127},{128,199}, // 189
{120,200},{201,115},{202,127},{128,203},{120,204},{205,115},{206,127}, // 196
{128,207},{120,208},{209,125},{210,127},{128,211},{130,212},{213,125}, // 203
{214,137},{138,215},{130,216},{217,125},{218,137},{138,219},{130,220}, // 210
{221,125},{222,137},{138,223},{130,224},{225,125},{226,137},{138,227}, // 217
{130,228},{229,125},{230,137},{138,231},{130,232},{233,125},{234,137}, // 224
{138,235},{130,236},{237,125},{238,137},{138,239},{130,240},{241,125}, // 231
{242,137},{138,243},{130,244},{245,135},{246,137},{138,247},{140,248}, // 238
{249,135},{250, 69},{ 80,251},{140,252},{249,135},{250, 69},{ 80,251}, // 245
{140,252},{  0,  0},{  0,  0},{  0,  0}};  // 252
#define nex(state,sel) State_table[state][sel]

这个表实际上给出了状态之间转换的方式。 比如说,初始状态是0,此时,如果下一个bit是1,那么就转换为State_table[0][1]这个状态,通过查表可知为2。 接着,如果再下一个bit是0,那么又转换为State_table[2][0]这个状态,通过查表可知为3。 所以,实际上当前的状态就表示了过去一定长度的bit流,但是这个长度是不定的。

LPAQ利用哈希表这样一个结构,保存了上一次出现某一个order 2, order 3, order 4, order 6, unigram的context时候,当时所处的状态。 这样,在我们处理新数据时,有如果再次遇到了同样的context,就获得上次出现这个context的状态,再把这个状态输入到对应的StateMap中(7种context每种都对应一个StateMap),以获得下一个bit是1的预测概率。

哈希表的代码如下所示

template<int B>
class HashTable {
    U8 *t;  // table: 1 element = B bytes: checksum priority data data
    const int N;  // size in bytes
public:
    HashTable(int n);
    
    U8 *operator[](U32 i);
    void clear();
    void copy(HashTable *bht) {
      memcpy(bht->t, t, sizeof(U8) * (N + B * 4));
    }
};
template<int B>
void HashTable<B>::clear() {
  memset(t, 0, sizeof(U8) * (N + B * 4));
  t += 64 - int(((long) t) & 63);
}
template<int B>
HashTable<B>::HashTable(int n): t(0), N(n) {
  assert(B >= 2 && (B & B - 1) == 0);
  assert(N >= B * 4 && (N & N - 1) == 0);
  alloc(t, N + B * 4 + 64);
  ori = t;
  t += 64 - int(((long) t) & 63);  // align on cache line boundary // cache line对齐
}
template<int B>
inline U8 *HashTable<B>::operator[](U32 i) {
  // 哈希计算,映射i到t上的某个位置。虽然实际上i本身就已经是order 1,2,3,4,6,unigram的context的哈希值。
  // 哈希表的entry为B个字节,使用哈希表的时候B为16,意思为entry为16字节。
  i *= 123456791;
  i = i << 16 | i >> 16; // 高低16bit互换
  i *= 234567891;
  int chk = i >> 24;
  i = i * B & N - B; // 把i映射到第i个16字节的位置。由于B为16,N-B正好使最后四个bit为0,保证最终计算出的i是16的整数倍。
  if (t[i] == chk) return t + i;  // 检查i的位置的checksum。如果相等,说明该context有历史记录。如果不相等,说明是哈希冲突或者从未写过这条记录。
  if (t[i ^ B] == chk) return t + (i ^ B); // 检查i^B位置的checksum。i^B相当于翻转倒数第五个bit(如果原来是1则变0,如果原来是0则变1)。i^B和i应该是相邻位置。
  if (t[i ^ B * 2] == chk)
    return t + (i ^ B * 2); // 检查i^B*2位置的checksum。i^B*2相当于翻转倒数第六个bit。i^B*2和i应该是隔了两个位置。 所以以上两个步骤相当于是哈希冲突时的线性探测法。
  // 如果全都冲突,或者没有记录,则尝试挑选出一个位置进行覆盖。
  if (t[i + 1] > t[i + 1 ^ B] || t[i + 1] > t[i + 1 ^ B * 2])
    i ^= B; // i+1的位置,代表i这个entry的state。根据计算顺序i+1^B实际上是(i+1)^B,表示对i+1的倒数第五个bit的翻转,相当于(i^B)+1,这里指向i^B位置entry的state。
  // 同样的道理,i+1^B*2指的是i^B*2这个位置的entry的state。
  // 这一行的含义是如果i位置的entry比i^B和i^B*2的entry的state都大,则把i置为i^B。这个前提应该是i的位置确实写有记录。
  if (t[i + 1] > t[i + 1 ^ B ^ B * 2]) i ^= B ^ B * 2;  // i+1^B^B*2表示对i同时反转倒数第五和倒数第六个bit。
  // 如果上一行运行成功,那么i+1^B^B*2就意味着把上一行的结果再翻转两次,相当于得到i+1^B*2。
  // 如果上一行运行失败,那么就会和i隔了三个位置。如果i处entry的state比i^B^B*2处entry的state还要大,则把i置为i^B^B*2
  // 这里的逻辑很奇怪。上一行的判断是带有或的,也就是说,如果第一个判断成功了,那么就直接给i赋值为(i^B),然后这一行再进行(i^B)和(i^B*2)之间的比较。
  // 如果上一行的第一个判断失败,第二个判断成功,那么这一行的判断就是注定失败的。
  // 如果上一行的判断全部失败,那么这一行的结果却不好说,有可能会将i置为(i^B^B*2),这会导致以后根本查询不到,这不就白写入了吗?
  memset(t + i, 0, B);   // 根据拟定的位置清除原有entry
  t[i] = chk; // 设置checksum
  return t + i; // 返回内容
}

// todo..