就地基数排序

Question

这是一个长文本。请多多包涵。归结起来，问题是：是否有可行的就地基数排序算法？

---

初步

我有大量的固定长度的字符串，它们只使用字母“A”、“C”、“G”和“T”（是的，你已经猜到了它：DNA) 我要排序。

目前，我使用std::sort，它在STL 的所有常见实现中使用introsort。这工作得很好。但是，我确信 radix sort 非常适合我的问题集，并且在实践中应该大大更好地工作。

详情

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

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

用例

理想情况下，此算法应适用于 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 字符串的就地基数排序？

Answer 1

虽然接受的答案完美地回答了问题的描述，但我已经到达这个地方，但徒劳无功地寻找一种将内联数组划分为 N 部分的算法。我自己写了一篇，就在这里。

警告：这不是一种稳定的分区算法，因此对于多级分区，必须重新分区每个结果分区而不是整个数组。优点是它是内联的。

它有助于解决所提出的问题的方式是，您可以根据字符串的字母重复内联分区，然后在分区足够小时使用您选择的算法对分区进行排序。

  function partitionInPlace(input, partitionFunction, numPartitions, startIndex=0, endIndex=-1) {
    if (endIndex===-1) endIndex=input.length;
    const starts = Array.from({ length: numPartitions + 1 }, () => 0);
    for (let i = startIndex; i < endIndex; i++) {
      const val = input[i];
      const partByte = partitionFunction(val);
      starts[partByte]++;
    }
    let prev = startIndex;
    for (let i = 0; i < numPartitions; i++) {
      const p = prev;
      prev += starts[i];
      starts[i] = p;
    }
    const indexes = [...starts];
    starts[numPartitions] = prev;
  
    let bucket = 0;
    while (bucket < numPartitions) {
      const start = starts[bucket];
      const end = starts[bucket + 1];
      if (end - start < 1) {
        bucket++;
        continue;
      }
      let index = indexes[bucket];
      if (index === end) {
        bucket++;
        continue;
      }
  
      let val = input[index];
      let destBucket = partitionFunction(val);
      if (destBucket === bucket) {
        indexes[bucket] = index + 1;
        continue;
      }
  
      let dest;
      do {
        dest = indexes[destBucket] - 1;
        let destVal;
        let destValBucket = destBucket;
        while (destValBucket === destBucket) {
          dest++;
          destVal = input[dest];
          destValBucket = partitionFunction(destVal);
        }
  
        input[dest] = val;
        indexes[destBucket] = dest + 1;
  
        val = destVal;
        destBucket = destValBucket;
      } while (dest !== index)
    }
    return starts;
  }

Answer 2

您当然可以通过以位为单位对序列进行编码来降低内存要求。 您正在查看排列，因此，对于长度为 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;
}

Answer 3

嗯，这是一个简单的 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 的方法，而不是通用的，但它应该很快。

编辑：

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

Answer 4

首先，考虑问题的编码。摆脱字符串，用二进制表示替换它们。使用第一个字节表示长度+编码。或者，在四字节边界处使用固定长度表示。然后基数排序变得容易得多。对于基数排序，最重要的是不要在内循环的热点进行异常处理。

好的，我对 4-nary 问题想得更多。为此，您需要像 Judy tree 这样的解决方案。下一个解决方案可以处理变长字符串；对于固定长度，只需删除长度位，这实际上更容易。

分配 16 个指针的块。可以重复使用指针的最低有效位，因为您的块将始终对齐。您可能需要一个特殊的存储分配器（将大存储分成更小的块）。有许多不同类型的块：

• 使用 7 位可变长度字符串进行编码。当它们填满时，您将它们替换为：
• 位置编码接下来的两个字符，您有 16 个指向下一个块的指针，以：
• 字符串最后三个字符的位图编码。

对于每种类型的块，您需要在 LSB 中存储不同的信息。由于您有可变长度的字符串，因此您也需要存储字符串结尾，并且最后一种块只能用于最长的字符串。随着您对结构的深入了解，应将 7 个长度位替换为 less。

这为您提供了一个相当快速且非常节省内存的排序字符串存储。它的行为有点像trie。要使其正常工作，请确保构建足够的单元测试。您想要覆盖所有块转换。您只想从第二种块开始。

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

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

您对处理双打感兴趣吗？对于短序列，将会有。调整块以处理计数很棘手，但它可以非常节省空间。

Answer 5

我会冒险并建议您切换到 heap/heapsort 实现。这个建议带有一些假设：

1. 您控制数据的读取
2. 只要“开始”排序，就可以对排序后的数据做一些有意义的事情。

heap/heap-sort 的美妙之处在于，您可以在读取数据的同时构建堆，并且可以在构建堆的那一刻开始获得结果。

让我们退后一步。如果你很幸运可以异步读取数据（即你可以发布某种读取请求并在某些数据准备好时收到通知），然后你可以在等待的同时构建一个堆下一块数据进入 - 甚至来自磁盘。通常，这种方法可以将一半排序的大部分成本隐藏在获取数据所花费的时间之后。

读取数据后，第一个元素已经可用。根据您发送数据的位置，这可能很棒。如果您将它发送到另一个异步阅读器，或者一些并行的“事件”模型或 UI，您可以随时发送块和块。

也就是说 - 如果您无法控制数据的读取方式，并且数据是同步读取的，并且在完全写出之前对已排序的数据没有用处 - 请忽略所有这些。 :(

查看维基百科文章：

• Heapsort
• Binary heap

Answer 6

“Radix sorting with no extra space”是一篇解决您的问题的论文。

Answer 7

你会想看看 Drs 的 Large-scale Genome Sequence Processing。笠原和森下。

由四个核苷酸字母 A、C、G 和 T 组成的字符串可以专门编码为整数，以便更快更快地处理。基数排序是本书讨论的众多算法之一。您应该能够调整这个问题的公认答案，并看到性能大幅提升。

Answer 8

在性能方面，您可能希望查看更通用的字符串比较排序算法。

目前你最终会触及每个字符串的每个元素，但你可以做得更好！

特别是，burst sort 非常适合这种情况。作为奖励，由于 Burstsort 基于尝试，它对于 DNA/RNA 中使用的小字母大小非常有效，因为您不需要在尝试实施。这些尝试也可能对您的类似后缀数组的最终目标有用。

在http://sourceforge.net/projects/burstsort/ 的源代码伪造上提供了一个不错的通用突发排序实现 - 但它不是就地的。

出于比较目的，http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf 中涵盖的 C-burstsort 实现在某些典型工作负载中的基准测试速度比快速排序和基数排序快 4-5 倍。

Answer 9

Radix-Sort 不考虑缓存，也不是对大集合最快的排序算法。 你可以看看：

• ti7qsort。 ti7qsort 是最快的整数排序（可用于固定大小的小字符串）。
• Inline QSORT
• String sorting

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

Answer 10

我会burstsort 字符串的压缩位表示。据称 Burstsort 比基数排序具有更好的局部性，使用突发尝试代替经典尝试来降低额外空间使用。原论文有尺寸。

Answer 11

看起来您已经解决了这个问题，但为了记录，似乎可行的就地基数排序的一个版本是“美国国旗排序”。它在这里描述：Engineering Radix Sort。一般的想法是对每个字符进行 2 次传递 - 首先计算每个字符有多少，因此您可以将输入数组细分为 bin。然后再次检查，将每个元素交换到正确的 bin 中。现在在下一个字符位置对每个 bin 进行递归排序。

Answer 12

dsimcha 的 MSB 基数排序看起来不错，但 Nils 更接近问题的核心，因为它观察到缓存局部性是在大问题规模时扼杀你的原因。

我建议一个非常简单的方法：

1. 根据经验估计基数排序有效的最大大小m。
2. 一次读取 m 元素块，对它们进行基数排序，然后将它们写出（如果有足够的内存，则写入内存缓冲区，否则写入文件），直到用完输入。
3. Mergesort 生成的排序块。

Mergesort 是我所知道的对缓存最友好的排序算法：“从数组 A 或 B 中读取下一项，然后将一项写入输出缓冲区。”它在磁带驱动器上高效运行。它确实需要2n 空间来对n 项目进行排序，但我敢打赌，您将看到大大改进的缓存位置将使这一点变得不重要——如果您使用的是非就地基数排序，您无论如何都需要额外的空间。

最后请注意，归并排序可以在没有递归的情况下实现，实际上这样做清楚了真正的线性内存访问模式。

Answer 13

您可以尝试使用trie。对数据进行排序只是遍历数据集并插入它；该结构是自然排序的，您可以将其视为类似于 B 树（除了不进行比较，您始终使用指针间接）。

缓存行为将有利于所有内部节点，因此您可能不会对此进行改进；但是您也可以摆弄 trie 的分支因子（确保每个节点都适合单个缓存行，分配类似于堆的 trie 节点，作为表示级别顺序遍历的连续数组）。由于尝试也是数字结构（长度为 k 的元素 O(k) 插入/查找/删除），因此您应该具有与基数排序相比具有竞争力的性能。

Answer 14

我从未见过就地基数排序，从基数排序的性质来看，只要临时数组适合内存，我怀疑它是否比非就地排序快得多。

原因：

排序对输入数组进行线性读取，但所有写入几乎都是随机的。从某个 N 向上，这归结为每次写入的缓存未命中。这种缓存未命中会减慢您的算法。无论是否到位都不会改变这种效果。

我知道这不会直接回答您的问题，但如果排序是一个瓶颈，您可能希望将 近似排序 算法作为 预处理步骤（软堆上的 wiki 页面可能会帮助您入门）。

这可以提供非常好的缓存局部性提升。教科书异地基数排序将表现更好。写入仍然几乎是随机的，但至少它们会聚集在相同的内存块周围，从而提高缓存命中率。

我不知道它在实践中是否可行。

顺便说一句：如果您只处理 DNA 字符串：您可以将一个 char 压缩成两位并大量打包您的数据。这将比简单的表示减少四倍的内存需求。寻址变得更加复杂，但无论如何，您的 CPU 的 ALU 在所有缓存未命中期间都需要花费大量时间。

Answer 15

如果您的数据集如此之大，那么我认为基于磁盘的缓冲区方法是最好的：

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 个存储桶），这样您就可以减少基于磁盘的缓冲区的分支。这可能会或可能不会提高性能，因此请尝试一下。