这是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万个序列,它比优化的内含子排序快3倍左右。

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

我知道这不会直接回答您的问题,但如果排序是一个瓶颈,您可能需要看看 近分选 算法作为 (软堆上的wiki页面可能会帮助您入门)。

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

顺便说一句:如果你只处理DNA字符串:你可以将一个字符压缩成两位,并大量打包数据。这将使一个虚拟表示的内存需求减少四倍。寻址变得更加复杂,但CPU的ALU在所有缓存未命中期间都有大量时间。

您当然可以通过以位编码序列来降低内存需求。 对于长度3,这是64个状态,可以用6位编码。所以序列中的每个字母看起来是2位,或者像你说的16个字符大约是32位。

因此,对于长度为3的序列,可以创建64个桶,大小可能为uint32或uint64。 将它们初始化为零。 将其用作下标,并递增该bucket。

按顺序遍历64个bucket,对于在该bucket中找到的计数,生成由该bucket表示的序列的多个实例。

一个4位的序列加上2位,因此将有256个存储桶。

在某个时刻,桶的数量将接近你的极限。

我认为这比在原地排序要快,因为桶很可能适合您的工作环境。

下面是一个展示这种技术的黑客

#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;
}

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

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的bucket(总共256个bucket),这样可以减少基于磁盘的缓冲区的分支。这可能会提高性能,也可能不会,所以请尝试一下。

我要冒险出去,建议你换成一堆/ heapsort 实施这一建议附带了一些假设:

1. 您可以控制数据的读取

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

让我们后退一步。如果您非常幸运,可以异步读取数据(也就是说,您可以发布某种读取请求,并在某些数据准备就绪时收到通知),那么您可以在等待下一个数据块进入时构建一个堆块,即使是从磁盘。通常,这种方法可以将排序的一半成本埋没在获取数据的时间之后。

请参阅维基百科文章:

• Heapsort
• Binary heap

" Radix sorting with no extra space

就性能而言,您可能希望了解更通用的字符串比较排序算法。

目前,你会接触到每个字符串的每个元素,但你可以做得更好!

burst sort

burstsort的一个体面的通用实现可在SourceForge上获得,网址为 http://sourceforge.net/projects/burstsort/ -但它还没有到位。

http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf 对于一些典型的工作负载,基准测试比快速排序和基数排序快4-5倍。

你会想看看 Large-scale Genome Sequence Processing

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

trie . 对数据进行排序只是对数据集进行迭代并插入它;结构是自然排序的,您可以将其视为类似于B-树(除非您不进行比较,而是 总是

缓存行为将有利于所有内部节点,因此您可能不会在这方面有所改进;但您也可以调整trie的分支因子(确保每个节点都适合单个缓存线,将类似于堆的trie节点分配为表示级别顺序遍历的连续数组)。由于try也是数字结构(长度为k的元素的O(k)insert/find/delete),因此与基数排序相比,您应该具有竞争性的性能。

我会的 burstsort

看起来你已经解决了这个问题,但是记录在案,一个可行的就地基数排序的版本是“美国国旗排序”。这里描述的是: Engineering Radix Sort . 一般的想法是对每个字符进行两次传递-首先计算每个字符的数量,以便将输入数组细分为多个存储单元。然后再次检查,将每个元素交换到正确的容器中。现在在下一个字符位置递归地对每个箱子排序。

你可以看看:

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

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

dsimcha的MSB基数排序看起来不错,但是Nils更接近问题的核心,因为它观察到缓存局部性在很大程度上是导致问题死亡的原因。

我建议采用一种非常简单的方法:

1. 根据经验估算最大尺寸 m
2. 读取数据块 M
3. 归并排序 生成的已排序块。

Mergesort是我所知道的对缓存最友好的排序算法:“从数组A或B中读取下一项,然后将一项写入输出缓冲区。”它在 磁带机 2n n

最后请注意,mergesort可以在没有递归的情况下实现,事实上,这样做可以清楚地说明真正的线性内存访问模式。

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

Judy tree 为了这个。下一个解决方案可以处理可变长度字符串;对于固定长度,只需删除长度位,这实际上使它更容易。

分配16个指针的块。指针的最低有效位可以重用,因为块总是对齐的。您可能需要一个特殊的存储分配器(将大型存储拆分为较小的块)。有许多不同类型的块:

• 使用可变长度字符串的7个长度位进行编码。当它们填满时,您可以将其替换为:
• 位置编码下两个字符,您有16个指向下一个块的指针,以:

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

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

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