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


初步

我有大量固定长度的小字符串,只使用字母“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字符串上的就地基数排序?


当前回答

虽然公认的答案完美地回答了问题的描述,但我已经到达了这个地方,徒劳地寻找一种算法将一个数组内联划分为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;
  }

其他回答

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

目前,您将触及每个字符串的每个元素,但您可以做得更好!

特别地,爆发排序非常适合这种情况。另外,由于burstsort是基于try的,它对于DNA/RNA中使用的小字母非常有效,因为你不需要在trie实现中构建任何形式的三元搜索节点、散列或其他三节点压缩方案。这些尝试对于类似后缀数组的最终目标也很有用。

burstsort的一个不错的通用实现可以在源代码forge (http://sourceforge.net/projects/burstsort/)上找到——但还不到位。

为了进行比较,http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf基准测试中涵盖的C-burstsort实现对于一些典型的工作负载比快速排序和基数排序快4-5倍。

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

原因:

排序对输入数组进行线性读取,但所有写入几乎都是随机的。从特定的N开始,这可以归结为每次写入的缓存丢失。这种缓存缺失会减慢你的算法。它是否到位并不会改变这种效果。

我知道这不能直接回答您的问题,但是如果排序是一个瓶颈,那么您可能想要看看作为预处理步骤的近似排序算法(软堆上的wiki-page可以让您开始)。

这可以很好地提高缓存的局部性。课本上的错位基数排序会表现得更好。写入仍然几乎是随机的,但至少它们会聚集在相同的内存块周围,这样就增加了缓存命中率。

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

顺便说一句:如果你只处理DNA字符串:你可以将一个字符压缩成两个比特,并打包大量数据。这将把内存需求减少到原始表示的四倍。寻址变得更加复杂,但无论如何,CPU的ALU在所有缓存丢失期间都有大量的时间。

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

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

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

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

您还可以使用压缩并将DNA的每个字母编码为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;
}