一尘不染

就地基数排序

algorithm

这是很长的文字。请多多包涵。归根结底，问题是： 是否存在可行的就地基数排序算法 ？

初步

我有很多 小的固定长度 字符串，它们只使用我要排序的字母“ A”，“ C”，“ G”和“
T”（是的，您猜对了：DNA）。

目前，我在STL的所有常见实现中都std::sort使用了introsort。这很好。但是，我坚信，基数排序非常适合我的问题，在实践中应该
会
更好。

细节

我已经用非常幼稚的实现测试了这个假设，对于相对较小的输入（大约10,000），这是正确的（至少要快两倍以上）。但是，当问题规模变大（ N >
5,000,000）时，运行时间将大大降低。

原因很明显：基数排序需要复制整个数据（实际上，在我的幼稚实现中不止一次）。这意味着我已经将〜4
GiB放入我的主内存中，这显然会降低性能。即使没有，我也负担不起这么大的内存，因为问题的大小实际上变得更大了。

用例

理想情况下，该算法应适用于2到100之间的任何字符串长度，适用于DNA以及DNA5（允许附加通配符“ N”），甚至适用于带有IUPAC

歧义码的 DNA
（导致16个不同的值）。但是，我意识到所有这些情况都无法解决，因此我对速度的提高感到满意。该代码可以动态决定要调度到哪个算法。

研究

不幸的是，维基百科上关于基数排序的文章是没有用的。关于就地变体的部分是完整的垃圾。在上基数NIST-
DADS部分排序旁边不存在的。有一篇听起来很有希望的论文，叫做“
高效自适应就地基数排序”，它描述了算法“
MSL”。不幸的是，这篇论文也令人失望。

特别是，有以下几点。

首先，该算法包含一些错误，并且有很多无法解释的地方。特别是，它没有详述递归调用（我只是假设它增加或减少了一些指针来计算当前的shift和mask值）。同时，它采用了功能dest_group并dest_address没有给出定义。我看不到如何有效地实现这些功能（也就是说，在O（1）中；至少dest_address是不平凡的）。

最后但并非最不重要的一点是，该算法通过将数组索引与输入数组中的元素交换来实现就位。显然，这仅适用于数值数组。我需要在字符串上使用它。当然，我可以拧紧强类型，然后继续假设内存可以容忍我存储不属于它的索引。但这仅在我可以将字符串压缩到32位内存（假设32位整数）的情况下有效。那只是16个字符（在16>
log（5,000,000）的那一刻，我们忽略它）。

其中一位作者的另一篇论文完全没有给出准确的描述，但它给出了MSL的运行时间为亚线性的情况，这是完全错误的。

回顾一下 ：是否有希望找到一个可行的参考实现，或者至少一个对DNA字符串起作用的就地基数排序的良好伪代码/描述？

阅读 281

2020-07-28

共1个答案

一尘不染

好吧，这是DNA的MSD基数排序的简单实现。它是用D语言编写的，因为这是我使用最多的语言，因此最不可能犯傻错误，但是可以轻松地将其翻译为其他语言。它就位，但是需要2 * seq.length通过数组。

void radixSort(string[] seqs, size_t base = 0) {
    if(seqs.length == 0)
        return;

    size_t TPos = seqs.length, APos = 0;
    size_t i = 0;
    while(i < TPos) {
        if(seqs[i][base] == 'A') {
             swap(seqs[i], seqs[APos++]);
             i++;
        }
        else if(seqs[i][base] == 'T') {
            swap(seqs[i], seqs[--TPos]);
        } else i++;
    }

    i = APos;
    size_t CPos = APos;
    while(i < TPos) {
        if(seqs[i][base] == 'C') {
            swap(seqs[i], seqs[CPos++]);
        }
        i++;
    }
    if(base < seqs[0].length - 1) {
        radixSort(seqs[0..APos], base + 1);
        radixSort(seqs[APos..CPos], base + 1);
        radixSort(seqs[CPos..TPos], base + 1);
        radixSort(seqs[TPos..seqs.length], base + 1);
   }
}

显然，这是特定于DNA的，而不是通用的，但是应该很快。

编辑：

我很好奇这个代码是否真的有效，所以我在等待自己的生物信息学代码运行时对其进行了测试/调试。上面的版本现在已经过实际测试并可以运行。对于1千万个5个碱基的序列，它比优化的introsort快3倍。

2020-07-28