使用后缀数组算法进行Burrows-Wheeler变换



我已经成功地为我正在编写的压缩测试台实现了BWT阶段(使用常规字符串排序)。我可以应用BWT,然后进行BWT逆变换,输出与输入匹配。现在我想加快使用后缀数组创建BW索引表的速度。我发现了两种相对简单的、据说很快的后缀数组创建O(n)算法,DC3和SA-IS,它们都带有C++/C源代码。我尝试使用源代码(在这里也可以找到开箱即用的编译SA-IS源代码),但未能获得正确的后缀数组/BWT索引表。以下是我所做的:

  1. T=输入数据,SA=输出后缀数组,n=T的大小,K=字母表大小,BWT=BWT索引表

  2. 我处理8位字节,但这两种算法都需要一个零字节形式的唯一sentinel/EOF标记(DC3需要3个,SA-IS需要一个),因此我将所有输入数据转换为32位整数,将所有符号增加1,并附加sentinel零字节。我是T.

  3. 我创建了一个整数输出数组SA(DC3的大小为n,KA-IS的大小为n+1)并应用算法。我得到的结果与排序BWT转换类似,但有些值是奇数(请参阅UPDATE 1)。此外,两种算法的结果略有不同。SA-IS算法在前面产生一个多余的索引值,因此所有结果都需要向左复制一个索引(SA[i]=SA[i+1])。

  4. 要将后缀数组转换为正确的BWT索引,我从后缀数组值中减去1,进行模运算,并应具有BWT索引(根据此):BWT[I]=(SA[I]-1)%n。

这是我提供SA算法并转换为BWT的代码。你应该能够或多或少地插入论文中的SA构造代码:

std::vector<int32_t> SuffixArray::generate(const std::vector<uint8_t> & data)
{
    std::vector<int32_t> SA;
    if (data.size() >= 2)
    {
        //copy data over. we need to append 3 zero bytes, 
        //as the algorithm expects T[n]=T[n+1]=T[n+2]=0
        //also increase the symbol value by 1, because the algorithm alphabet is [1,K]
        //(0 is used as an EOF marker)
        std::vector<int32_t> T(data.size() + 3, 0);
        std::copy(data.cbegin(), data.cend(), T.begin());
        std::for_each(T.begin(), std::prev(T.end(), 3), [](int32_t & n){ n++; });
        SA.resize(data.size());
        SA_DC3(T.data(), SA.data(), data.size(), 256);
        OR
        //copy data over. we need to append a zero byte, 
        //as the algorithm expects T[n-1]=0 (where n is the size of its input data)
        //also increase the symbol value by 1, because the algorithm alphabet is [1,K] 
        //(0 is used as an EOF marker)
        std::vector<int32_t> T(data.size() + 1, 0);
        std::copy(data.cbegin(), data.cend(), T.begin());
        std::for_each(T.begin(), std::prev(T.end(), 1), [](int32_t & n){ n++; });
        SA.resize(data.size() + 1); //crashes if not one extra byte at the end
        SA_IS((unsigned char *)T.data(), SA.data(), data.size() + 1, 256, 4); //algorithm expects size including sentinel
        std::rotate(SA.begin(), std::next(SA.begin()), SA.end()); //rotate left by one to get same result as DC3
        SA.resize(data.size());
    }
    else
    {
        SA.push_back(0);
    }
    return SA;
}
void SuffixArray::toBWT(std::vector<int32_t> & SA)
{
    std::for_each(SA.begin(), SA.end(), [SA](int32_t & n){ n = ((n - 1) < 0) ? (n + SA.size() - 1) : (n - 1); });
}

我做错了什么?

更新1
当将算法应用于少量的测试文本数据时,如"yabbadabbado"/"这是一个测试。"/"abaaba"或大文本文件(坎特伯雷语料库中的alice29.txt),它们运行良好。实际上,toBWT()函数甚至没有必要
当将算法应用于包含完整8位字节字母表(可执行文件等)的文件中的二进制数据时,它们似乎无法正常工作。将算法的结果与常规BWT索引的结果进行比较,我注意到前面有错误的索引(在我的情况下是4)。索引的数量(偶然?)对应于算法的递归深度。索引指向原始源数据最后出现0的位置(在构建T时我将其转换为1之前)。。。

更新2
当我二进制比较常规BWT数组和后缀数组时,会有更多不同的值。这可能是意料之中的事,因为afair排序不一定与标准排序相同,但数组转换的结果数据应该相同。事实并非如此。

更新3
我尝试修改一个简单的输入字符串,直到两个算法都"失败"。将字符串"this is a test."的两个字节更改为255或0(从74686973206973206120746573742Eh更改为例如7468697320069732061FF7.465774FFh,最后一个字节必须更改!)索引和转换后的字符串不再正确。将字符串的最后一个字符更改为字符串中已经存在的字符似乎也足够了,例如"这是一个测试"746869732069732061207465737473h。然后,转换字符串的两个索引和两个字符将被交换(比较常规排序BWT和使用SA的BWT)。

我发现必须将数据转换为32位的整个过程有点尴尬。如果有人有更好的解决方案(纸上的,更好的是,一些源代码),可以直接从256个字符的字符串中生成后缀数组,我会很高兴。

我现在已经弄清楚了。我的解决方案有两个方面。有些人建议使用一个图书馆,我做了森裕太的SAIS lite
真正的解决方案是复制并连接输入字符串,并在此字符串上运行SA生成。保存输出字符串时,需要过滤掉所有超过原始数据大小的SA索引。这不是一个理想的解决方案,因为您需要分配两倍的内存,复制两次,并对两倍的数据量进行转换,但它仍然比std::sort快50-70%。如果你有更好的解决方案,我很乐意听到。
你可以在这里找到更新的代码。

最新更新