205

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

初步的

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

目前，我使用which 在STL的所有常见实现中std::sort使用introsort。这工作得很好。但是，我相信基数排序非常适合我的问题集，并且在实践中应该工作得更好。

细节

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

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

用例

理想情况下，该算法应适用于 2 到 100 之间的任何字符串长度，适用于 DNA 和 DNA5（允许额外的通配符“N”），甚至是具有IUPAC 歧义码的 DNA （产生 16 个不同的值）。但是，我意识到无法涵盖所有​​这些情况，因此我对获得的任何速度改进感到满意。代码可以动态决定分派到哪个算法。

研究

不幸的是，关于基数排序的维基百科文章毫无用处。关于就地变体的部分完全是垃圾。NIST-DADS 关于基数排序的部分几乎不存在。有一篇听起来很有前途的论文叫做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 字符串的就地基数排序？

4

15 回答 15

62

好吧，这是 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 倍。

于 2009-01-20T21:19:25.870 回答
21

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

原因：

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

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

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

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

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

于 2009-01-20T21:36:30.173 回答
9

您当然可以通过以位编码序列来降低内存要求。您正在查看排列，因此，对于长度为 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;
}
于 2009-01-20T22:48:02.010 回答
6

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

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

于 2009-01-20T21:24:34.970 回答
6

我将冒险出去并建议您切换到堆/堆排序实现。这个建议带有一些假设：

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

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

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

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

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

请参阅维基百科文章：

于 2009-01-20T22:00:44.570 回答
5

“没有多余空间的基数排序”是一篇解决您的问题的论文。

于 2010-08-05T07:07:06.063 回答
4

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

目前你最终会触及每根弦的每一个元素，但你可以做得更好！

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

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

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

于 2009-07-15T00:37:54.477 回答
4

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

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

于 2010-01-23T18:17:44.187 回答
3

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

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

于 2009-01-21T05:05:16.630 回答
3

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

于 2009-01-24T22:11:30.693 回答
2

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

于 2009-01-23T23:50:35.627 回答
2

Radix-Sort 没有缓存意识，也不是大集合最快的排序算法。你可以看看：

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

于 2009-06-14T10:25:39.277 回答
1

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

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

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

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

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

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

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

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

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

于 2009-01-20T22:45:21.907 回答
1

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

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

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

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

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

于 2009-01-21T11:40:03.013 回答
0

虽然接受的答案完美地回答了问题的描述，但我已经到了这个地方，徒劳地寻找一种将内联数组划分为 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;
  }
于 2020-09-06T10:48:04.980 回答