这是一篇很长的文章。请原谅我。归结起来,问题是:是否存在可行的就地基数排序算法?


初步

我有大量固定长度的小字符串,只使用字母“a”,“C”,“G”和“T”(是的,你已经猜到了:DNA),我想对它们进行排序。

目前,我使用std::sort,它在STL的所有常见实现中使用introsort。这工作得很好。然而,我确信基数排序完全适合我的问题集,在实践中应该工作得更好。

细节

我用一个非常简单的实现测试了这个假设,对于相对较小的输入(大约10,000),这是正确的(至少快两倍多)。然而,当问题规模变大(N > 5,000,000)时,运行时间会急剧下降。

原因很明显:基数排序需要复制整个数据(实际上在我的简单实现中不止一次)。这意味着我在主存中放置了~ 4 GiB,这显然会降低性能。即使它没有,我也不能使用这么多内存,因为问题的大小实际上会变得更大。

用例

理想情况下,该算法应该适用于2到100之间的任何字符串长度,适用于DNA和DNA5(允许额外的通配符“N”),甚至适用于具有IUPAC歧义码的DNA(导致16个不同的值)。然而,我意识到所有这些情况都无法涵盖,所以我对我得到的任何速度改进都很满意。代码可以动态地决定向哪个算法分派。

研究

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

具体来说,有以下几点。

首先,该算法包含了一些错误,并留下了许多无法解释的地方。特别是,它没有详细说明递归调用(我只是假设它增加或减少一些指针来计算当前的移位和掩码值)。此外,它使用函数dest_group和dest_address,但没有给出定义。我不知道如何有效地实现这些(也就是说,在O(1);至少dest_address不是简单的)。

Last but not least, the algorithm achieves in-place-ness by swapping array indices with elements inside the input array. This obviously only works on numerical arrays. I need to use it on strings. Of course, I could just screw strong typing and go ahead assuming that the memory will tolerate my storing an index where it doesn’t belong. But this only works as long as I can squeeze my strings into 32 bits of memory (assuming 32 bit integers). That's only 16 characters (let's ignore for the moment that 16 > log(5,000,000)).

另一篇论文的作者没有给出准确的描述,但它给出了MSL的运行时是次线性的,这是完全错误的。

回顾一下:有没有希望找到一个工作的参考实现,或者至少是一个好的伪代码/描述,一个工作在DNA字符串上的就地基数排序?


当前回答

基数排序不是缓存意识,也不是对大集最快的排序算法。 你可以看看:

ti7qsort。Ti7qsort是对整数最快的排序(可用于固定大小的小字符串)。 内联QSORT 字符串排序

您还可以使用压缩并将DNA的每个字母编码为2位,然后存储到排序数组中。

其他回答

这是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的,而不是普遍的,但它应该很快。

编辑:

我很好奇这个代码是否真的有效,所以在等待我自己的生物信息学代码运行时,我对它进行了测试/调试。上面的版本现在已经经过测试并且可以工作。对于1000万个每个5个碱基的序列,它比优化的引入排序快大约3倍。

首先,考虑问题的编码。去掉字符串,用二进制表示代替它们。使用第一个字节表示长度+编码。或者,在四字节边界上使用固定长度的表示。基数排序就简单多了。对于基数排序,最重要的是不要在内部循环的热点处进行异常处理。

好的,我想了一下4元的问题。你想要一个像朱迪树一样的解决方案。下一个解决方案可以处理可变长度的字符串;对于固定长度,只要去掉长度位,这实际上让它更简单。

分配16个指针的块。指针中最不重要的部分可以被重用,因为你的块总是对齐的。您可能需要为它使用一个特殊的存储分配器(将大的存储分解为较小的块)。有许多不同类型的积木:

用可变长度字符串的7个长度位进行编码。当它们填满时,你用: 位置编码接下来的两个字符,你有16个指针指向下一个块,以: 字符串最后三个字符的位图编码。

对于每种类型的块,您需要在lsb中存储不同的信息。当你有可变长度的字符串时,你也需要存储字符串的结尾,最后一种块只能用于最长的字符串。随着结构的深入,长度为7的位应该被更少的位所取代。

这为您提供了一个合理快速和非常有效的内存存储排序字符串。它会表现得有点像一个trie。要让它工作,请确保构建足够的单元测试。您希望覆盖所有块转换。你只想从第二种积木开始。

为了获得更好的性能,您可能需要添加不同的块类型和更大的块大小。如果块总是相同的大小和足够大,您可以为指针使用更少的位。块大小为16个指针,在32位地址空间中已经有一个字节空闲。查看Judy树文档,了解有趣的块类型。基本上,您添加代码和工程时间以进行空间(和运行时)权衡

您可能希望从头四个字符的256宽直接基数开始。这提供了一个不错的空间/时间权衡。在这个实现中,你得到的内存开销比简单的trie少得多;它大约小了三倍(我还没有测量过)。如果常数足够低,O(n)不是问题,就像你在与O(n log n)快速排序比较时注意到的那样。

你对处理双数感兴趣吗?对于短序列,会有。调整块来处理计数是很棘手的,但它可以非常节省空间。

你会想看一看博士们的大规模基因组序列处理。笠原和森下。

由四个核苷酸字母A、C、G和T组成的字符串可以特别编码为整数,以便更快地处理。基数排序是书中讨论的许多算法之一;您应该能够适应这个问题的公认答案,并看到一个很大的性能改进。

看起来您已经解决了这个问题,但是为了记录,似乎有一种可行的就地基数排序是“美国国旗排序”。这里描述了:工程基数排序。一般的想法是对每个字符进行2次传递——首先计算每个字符有多少个,这样就可以将输入数组细分为箱子。然后再执行一遍,将每个元素交换到正确的bin中。现在递归地对每个箱子的下一个字符位置排序。

你可以试着用trie。对数据排序就是简单地遍历数据集并插入数据集;结构是自然排序的,你可以把它想象成类似于b -树(除了不做比较,你总是使用指针指向)。

缓存行为将有利于所有内部节点,因此您可能无法在此基础上进行改进;但是您也可以修改树的分支因子(确保每个节点都适合单个缓存线,将树节点分配为类似于堆的节点,作为表示级别顺序遍历的连续数组)。由于try也是数字结构(对于长度为k的元素,O(k)插入/查找/删除),因此您的性能应该优于基数排序。