原地基数排序

Question

这是一段很长的文本，请耐心阅读。简而言之，问题是：是否有可行的原地基数排序算法？

---

初步

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

目前，我使用std::sort，它在所有常见的STL实现中使用introsort。这个方法效果还不错。然而，我相信基数排序完全适合我的问题集，并且在实践中应该表现得更好。

细节

我用一个非常天真的实现来测试了这个假设，对于相对较小的输入（大约为10,000），这是正确的（至少快两倍）。然而，当问题规模变得更大（N > 5,000,000）时，运行时间会极度下降。

原因很明显：基数排序需要复制整个数据（实际上在我的天真实现中不止一次）。这意味着我已经将~4 GiB放入了我的主内存中，这显然会影响性能。即使没有这个问题，我也无法承受使用这么多内存，因为问题的规模实际上会更大。

用例

Ideally，这个算法应该适用于2到100个字符长度的DNA和DNA5（允许使用额外的通配符“N”），甚至包括带有IUPAC ambiguity codes的DNA（导致16个不同的值）。然而，我意识到并不是所有情况都能覆盖，所以我很高兴能获得任何速度上的提升。代码可以动态地决定派发哪种算法。

研究

不幸的是，Wikipedia article on radix sort对原地变体的部分完全没有用处。NIST-DADS section on radix sort几乎不存在。有一篇听起来很有前途的论文，叫做Efficient Adaptive In-Place Radix Sorting，描述了算法“MSL”。不幸的是，这篇论文也令人失望。

特别是以下几点。

首先，该算法存在几个错误并留下了很多未解释的问题。特别是，它没有详细说明递归调用（我只是假设它增加或减少某些指针来计算当前的移位和掩码值）。此外，它在不给出定义的情况下使用函数dest_group和dest_address。我无法看出如何高效地实现这些（也就是说，在O（1）中；至少dest_address不是微不足道的）。

最后，该算法通过交换数组索引与输入数组内部的元素来实现原地性。这显然仅适用于数字数组。我需要在字符串上使用它。当然，我可以放弃强类型并继续假设内存将容忍我在不属于它的位置存储索引。但是，只要我能将我的字符串压缩到32位内存中（假设32位整数），这种方法才有效。那只有16个字符（暂且忽略16>log（5,000,000））。

该作者的另一篇论文根本没有提供准确的描述，但它将MSL的运行时间描述为次线性，这是完全错误的。

总之：是否有希望找到一个适用于DNA字符串的工作参考实现或至少是工作在原地的基数排序的好伪代码/描述？

- Konrad Rudolph

---

这是一个写得非常好的问题。 - JustinT

小固定长度字符串有多小？ - EvilTeach

@Stephan：这一切都很好。但是在复制/缓存未命中的情况下，我只会遇到延迟。而在内存方面，我会遇到物理限制。这是绝对不可谈判的。所有那些将数据的部分存储在磁盘上的花哨技术肯定比当前的快速排序解决方案慢。 - Konrad Rudolph

dsimcha的解决方案则在某些输入情况下绝对比quicksort更快。移动次数可能会很高，缓存局部性也较小，但在现实世界中仍然表现良好。我还稍微调整了一下解决方案，以减少需要执行的交换次数。 - Konrad Rudolph

@PeterMortensen 对于未来，我欣赏您提供更正和链接以增加上下文。然而，我并不特别欣赏那些仅仅是风格问题的编辑（如“在...次序中”与“按...顺序”，“即”与“也就是”）。 - Konrad Rudolph

Answer 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的具体内容，而不是一般性的内容，但速度应该很快。

编辑： 我很好奇这段代码是否有效，所以在等待我的生物信息学代码运行时进行了测试/调试。上面的版本现在已经经过测试并且能够正常工作。对于5000万个长度为5的序列，它比优化后的introsort快约3倍。

- dsimcha

---

如果您可以接受两次遍历的方法，这种方法也可以扩展到基数为N的情况：第一次遍历只需计算每个数字出现的次数。如果您正在对数组进行分区，则此操作会告诉您每个数字在哪里开始。第二次遍历则将该数字交换到数组中适当的位置。 - Jason S

这个基数排序似乎是美国国旗排序的一个特例，后者是一种广为人知的原地基数排序变体。 - Edward Kmett

完全没有关系，但是你用什么IDE来写D语言？一直在寻找一个好的IDE... - mpen

CodeBlocks很糟糕（基本上只是一个编辑器加上基本的构建自动化），我一直在寻找更好的替代品，但目前还没有适用于D2（该语言的最新版本）的好的集成开发环境。D语言具有足够的语言级别特性来消除不必要的样板代码，因此相比其他语言，你并不需要像其他语言那样需要一个好的IDE。 - dsimcha

这段代码更类似于具有预定义枢轴的快速排序。Jason S.的建议似乎更容易扩展到更长的前缀。 - qsantos

Answer 2

我从未见过原地基数排序，在基数排序的性质方面，只要临时数组可以放入内存中，我怀疑它比就地排序快不了多少。

原因：

排序对输入数组进行线性读取，但所有写操作几乎是随机的。从某个N开始，这归结为每次写入一个缓存未命中。这个缓存未命中会降低算法的速度。无论是否就地排序都不会改变这种影响。

我知道这不会直接回答你的问题，但如果排序是瓶颈，你可能需要查看作为预处理步骤的近似排序算法（维基百科上的软堆页面可能会为您提供帮助）。

这可以大大增加缓存本地性。一本教科书式的就地基数排序将表现更好。写操作仍然几乎是随机的，但至少它们将聚集在相同的内存块周围，并因此增加缓存命中率。

我不知道它在实践中是否有效。

顺便说一句：如果只处理DNA字符串：您可以将一个字符压缩为两位并相当紧凑地打包数据。这将使内存需求降低四分之一，远高于幼稚的表示法。地址寻址变得更加复杂，但是你的CPU ALU在所有缓存未命中期间有大量时间可以使用。

- Nils

---

两个好的观点；近排序是个新概念，我得去了解一下。缓存未命中是另一个需要考虑的问题，让我有些担心。;-) 我要看看这个。 - Konrad Rudolph

这对我来说也是新的（几个月），但一旦你掌握了概念，就会开始看到性能改进的机会。 - Nils Pipenbrinck

除非你的基数非常大，否则这些写入远非“几乎随机”。例如，假设你一次只对一个字符进行排序（基数为4），所有的写入都将被写入4个线性增长的桶之一。这既可以缓存，也可以预取。当然，你可能想要使用更大的基数，在缓存和预取友好性与基数大小之间存在权衡。你可以通过软件预取或为桶分配一个临时区域并定期刷新到“真实”桶中来将均衡点推向更大的基数。 - BeeOnRope

Answer 3

您可以通过将序列编码为位来降低内存要求。对于长度为2的排列，使用“ACGT”，有16个状态或4位。对于长度为3的排列，有64个状态，可以用6位编码。因此，每个字母序列需要2位，像您所说的那样，16个字符需要大约32位。

如果有一种方法可以减少有效单词的数量，可能还可以进行进一步压缩。

因此，对于长度为3的序列，可以创建64个桶，大小为uint32或uint64。将它们初始化为零。迭代通过您非常庞大的3个字符序列列表，并按上述方式对它们进行编码。使用这个作为下标并递增该桶。重复此过程，直到处理完所有序列。

接下来，重新生成您的列表。

按顺序迭代64个桶，对于在该桶中找到的计数，生成表示该桶的序列的实例。
当所有桶都被迭代时，您就有了排序后的数组。

长度为4的序列增加了2位，因此将有256个桶。长度为5的序列增加了2位，因此将有1024个桶。

在某些情况下，桶的数量将接近您的限制。如果从文件中读取序列而不是将它们保留在内存中，则可用于桶的内存将更多。

我认为这比在原地进行排序要快，因为桶很可能适合您的工作集内。以下是演示该技术的hack。

#include <iostream>
#include <iomanip>

#include <math.h>

using namespace std;

const int width = 3;
const int bucketCount = exp(width * log(4)) + 1;
      int *bucket = NULL;

const char charMap[4] = {'A', 'C', 'G', 'T'};

void setup
(
    void
)
{
    bucket = new int[bucketCount];
    memset(bucket, '\0', bucketCount * sizeof(bucket[0]));
}

void teardown
(
    void
)
{
    delete[] bucket;
}

void show
(
    int encoded
)
{
    int z;
    int y;
    int j;
    for (z = width - 1; z >= 0; z--)
    {
        int n = 1;
        for (y = 0; y < z; y++)
            n *= 4;

        j = encoded % n;
        encoded -= j;
        encoded /= n;
        cout << charMap[encoded];
        encoded = j;
    }

    cout << endl;
}

int main(void)
{
    // Sort this sequence
    const char *testSequence = "CAGCCCAAAGGGTTTAGACTTGGTGCGCAGCAGTTAAGATTGTTT";

    size_t testSequenceLength = strlen(testSequence);

    setup();


    // load the sequences into the buckets
    size_t z;
    for (z = 0; z < testSequenceLength; z += width)
    {
        int encoding = 0;

        size_t y;
        for (y = 0; y < width; y++)
        {
            encoding *= 4;

            switch (*(testSequence + z + y))
            {
                case 'A' : encoding += 0; break;
                case 'C' : encoding += 1; break;
                case 'G' : encoding += 2; break;
                case 'T' : encoding += 3; break;
                default  : abort();
            };
        }

        bucket[encoding]++;
    }

    /* show the sorted sequences */ 
    for (z = 0; z < bucketCount; z++)
    {
        while (bucket[z] > 0)
        {
            show(z);
            bucket[z]--;
        }
    }

    teardown();

    return 0;
}

- EvilTeach

---

为什么要比较，当你可以哈希呢？ - wowest

当然没错。 性能通常是任何DNA处理中的一个问题。 - EvilTeach

Answer 4

我建议你切换到堆/堆排序实现。这个建议带有一些假设：

1. 你控制数据的读取
2. 你可以在开始获取排序后的数据时立即对其进行有意义的操作。

堆/堆排序的美妙之处在于，你可以在读取数据的同时构建堆，而且一旦构建好堆就可以立即开始获取结果。

让我们退后一步。如果你很幸运可以异步地读取数据（也就是说，你可以发布某种类型的读取请求，并在准备好一些数据时得到通知），那么你可以在等待下一批数据到来的同时构建一部分堆——即使是从磁盘中读取。通常情况下，这种方法可以将排序的大部分成本隐藏在获取数据的时间背后。

一旦你读取了数据，第一个元素就已经可用了。根据你发送数据的位置，这可能非常好。如果你将其发送到另一个异步读取器、某个并行的“事件”模型或者 UI 中，你可以随着进展不断地发送块和块的数据。

话虽如此 - 如果您无法控制数据的读取方式，并且它是同步读取的，并且在完全写出排序数据之前您没有使用排序数据的用途 - 那么请忽略所有这些。 :(

参见维基百科文章：

• 堆排序
• 二叉堆

- Joe

---

好的建议。然而，在我的特定情况下，维护堆的开销比仅在向量中累积数据并在所有数据到达后进行排序要大。 - Konrad Rudolph

Answer 5

如果您的数据集非常大，那么我认为最好采用基于磁盘缓存的方法：

sort(List<string> elements, int prefix)
    if (elements.Count < THRESHOLD)
         return InMemoryRadixSort(elements, prefix)
    else
         return DiskBackedRadixSort(elements, prefix)

DiskBackedRadixSort(elements, prefix)
    DiskBackedBuffer<string>[] buckets
    foreach (element in elements)
        buckets[element.MSB(prefix)].Add(element);

    List<string> ret
    foreach (bucket in buckets)
        ret.Add(sort(bucket, prefix + 1))

    return ret

我建议你尝试将数据分成更多的桶，例如，如果你的字符串是：

GATTACA

第一个 MSB 调用会返回 GATT 的桶（共 256 个桶），这样可以减少基于磁盘缓冲区的分支。这可能会提高性能，因此请进行实验。

- FryGuy

---

我们在某些应用程序中使用内存映射文件。然而，一般情况下，我们的工作基于这样的假设：机器提供的RAM仅仅足够不需要显式磁盘备份（当然，交换仍然会发生）。但我们已经在开发自动磁盘备份数组的机制。 - Konrad Rudolph

Answer 6

"无需额外空间的基数排序" 是一篇解决你问题的论文。"

- eig

---

看起来很有前途，虽然这个问题已经被解决了。不过，这会被加入我的参考库。 - Konrad Rudolph

Answer 7

就性能而言，您可能需要查看更通用的字符串比较排序算法。

目前，您会接触到每个字符串的每个元素，但您可以做得更好！

特别是，burst sort非常适合这种情况。作为奖励，由于burstsort基于tries，因此对于DNA / RNA中使用的小字母表大小非常有效，因为您不需要将任何三进制搜索节点，哈希或其他trie节点压缩方案构建到trie实现中。这些tries也可以用于您的后缀数组式最终目标。

在source forge上提供了一个不错的burstsort通用实现，网址为http://sourceforge.net/projects/burstsort/，但它不是原地排序。

为了比较，C-burstsort实现在http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf中涵盖的基准测试比快速排序和基数排序快4-5倍。

- Edward Kmett

---

我肯定要研究一下Burst Sort - 尽管目前我不知道如何建立原位Trie。在实际应用中，后缀数组已经几乎取代了后缀树（因此也就是Trie）在生物信息学中，因为其性能优越。 - Konrad Rudolph

Answer 8

您需要查看Kasahara博士和Morishita博士的大规模基因组序列处理。该书中讨论了许多算法，包括可以将由四个核苷酸字母A、C、G和T组成的字符串特殊编码为整数，以实现更快速的处理。其中之一是基数排序，您可以适应此问题的已接受答案，并看到性能大幅提高。

- Rudiger

---

本书介绍的基数排序不是就地排序，因此不能用于此目的。至于字符串压缩，我当然已经做到了这一点。我的最终解决方案（以下发布）没有显示这一点，因为库使我可以像处理普通字符串一样处理它们-但是使用的“RADIX”值当然可以调整为更大的值。 - Konrad Rudolph

Answer 9

您可以尝试使用Trie。数据排序只需遍历数据集并插入即可；该结构自然排序，可以将其视为类似于B树（除了进行比较外，您始终使用指针间接引用）。缓存行为将支持所有内部节点，因此您可能无法改进它；但是您也可以调整trie的分支因子（确保每个节点适合单个缓存行，分配trie节点类似于堆，作为代表级别顺序遍历的连续数组）。由于tries也是数字结构（O(k)插入/查找/删除长度为k的元素），因此您应该具有与基数排序相竞争的性能。

- Tom

---

字典树与我天真的实现一样，存在相同的问题：它需要额外的O(n)内存，这太多了。 - Konrad Rudolph

Answer 10

我会使用burstsort对字符串进行压缩位表示。据称，Burstsort比基数排序具有更好的局部性，在经典Tries的位置上使用Burst Tries可以将额外的空间使用降至最低。原始论文中有测量数据。

- Darius Bacon