就地基数sorting
这是一个长文本。 请多多包涵。 下来，问题是： 是否有一个可行的就地基数sortingalgorithm ？
初步
我有很多小的固定长度的string，只使用我想分类的字母“A”，“C”，“G”和“T”（是的，你猜对了）。
目前，我使用std::sort ，它在STL的所有常见实现中使用introsort 。 这工作得很好。 不过，我确信基数sorting完全适合我的问题，在实践中应该更好。
细节
我已经用非常幼稚的实现来testing这个假设，对于相对较小的input（大约为10,000），这是真实的（至less是两倍以上）。 但是，当问题规模变大（ N > 5,000,000）时，运行时间会严重恶化。
原因很明显：基数sorting需要复制整个数据（实际上我的幼稚实现不止一次）。 这意味着我已经将〜4 GiB放入了我的主内存，这显然会导致性能下降。 即使没有，我也不能使用这么多的记忆，因为问题的规模实际上变得更大了。
用例
理想的情况下，这个algorithm应该适用于DNA和DNA5（它允许一个额外的通配符“N”），甚至是带有IUPAC 模糊代码的 DNA（导致16个不同的值）。 但是，我意识到所有这些情况都不能被覆盖，所以我对任何速度的改善感到满意。 代码可以dynamic地决定分派给哪个algorithm。
研究
不幸的是， 关于基数sorting的维基百科文章是无用的。 关于就地变体的部分是完整的垃圾。 基数sorting的NIST-DADS部分几乎不存在。 有一种看起来很有前途的论文，称为高效自适应就地基sorting ，它描述了“MSL”algorithm。 不幸的是，这篇论文也令人失望。
具体来说，有以下几点。
首先，该algorithm包含几个错误，并留下很多原因。 特别是，它没有详细说明recursion调用（我简单地假设它递增或减less一些指针来计算当前的移位和掩码值）。 而且，它使用函数dest_group和dest_address而不给定义。 我没有看到如何有效地实现这些（即在O（1）;至lessdest_address是不平凡的）。
最后但并非最不重要的是，该algorithm通过将数组索引与input数组内的元素交换来实现就地性。 这显然只适用于数值数组。 我需要在string上使用它。 当然，我可以拧紧打字，然后继续前进，假设内存将容忍我存储一个不属于它的索引。 但是，只要我能把我的string压缩到32位内存（假设是32位整数），这个工作就行得通。 那只有16个字符（让我们暂且忽略16> log（5,000,000））。
其中一位作者的另一篇论文根本没有给出准确的描述，但是它使得MSL的运行时呈现出平坦的错误。
回顾一下 ：有没有希望find一个工作的参考实现，或至less有一个良好的伪代码/描述工作在地方基数sorting的DNAstring？
那么，这里是一个简单的DNA的MSD基数sorting实现。 它是用D编写的，因为这是我使用最多的语言，因此最不可能犯的错误，但它可以很容易地翻译成其他语言。 它在原地，但需要2 * seq.length穿过arrays。
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的具体情况，而不是一般的，但它应该是快速的。
编辑：我很好奇这个代码是否真正起作用，所以我在等待我自己的生物信息学代码运行时对其进行了testing/debugging。 现在上面的版本实际上是testing和工作。 对于每个5个碱基的1000万个序列，它比优化的插入片段快大约3倍。
我从来没有见过一种就地的基数sorting，从基数sorting的性质来看，只要临时数组适合内存，我怀疑它比不合适的sorting要快得多。
原因：
sorting在input数组上进行线性读取，但是所有的写入操作几乎是随机的。 从某个N向上，这可以归结为每个写入的caching缺失。 这个caching未命中是减慢你的algorithm。 如果有的话不会改变这个效果。
我知道这不会直接回答你的问题，但如果sorting是一个瓶颈，你可能需要看看sortingalgorithm作为预处理步骤 （软堆上的维基页面可能让你开始）。
这可以给一个非常好的caching局部性提升。 然后，一个不符合地方的文本书本将会执行得更好。 写入仍然是随机的，但是至less它们将围绕相同的内存块聚集，并且因此增加caching命中率。
我不知道它是否在实践中成功。
顺便说一句：如果你只处理DNAstring：你可以压缩一个字符两位，并打包你的数据相当多。 这将减less内存需求的四分之一以上的代表性。 寻址变得越来越复杂，但是CPU的ALU在所有的caching未命中期间都有大量的时间花费。
基于dsimcha的代码，我实现了一个更适合我们使用的框架（SeqAn）的更通用的版本。 实际上，移植代码是非常简单的。 后来才发现有关这个话题的刊物。 伟大的事情是：他们基本上和你们一样说话。 Andersson和Nilsson关于实现Radixsort的论文绝对值得一读。 如果你碰巧知道德语，那么一定要读David Weese的gradle论文 ，他会实现一个通用的子串索引。 本文的大部分内容都是致力于详细分析构build索引的成本，考虑二级存储和极大文件。 他的工作成果实际上是在SeqAn实施的，而不是在我需要的地方。
只是为了好玩，这里是我写的代码（我不认为有人不使用SeqAn会有任何用处）。 注意，它仍然不考虑基数大于4.我预计这会对性能产生巨大的影响，但不幸的是我现在根本没有时间去实现这一点。
代码的执行速度是短string的Introsort的两倍以上。 盈亏平衡点在12-13左右。 string的types（例如，是否具有4,5或16个不同的值）是相对不重要的。 从人类基因组的第2号染色体上分离> 6,000,000个DNA读取在我的PC上花费超过2秒钟。 只是为了logging，这很快 ！ 特别是考虑到我不使用SIMD或任何其他硬件加速。 此外，valgrind指出我的主要瓶颈是operator new在string赋值中是operator new的。 它被称为约65,000,000次 – 每个string十次！ 对于这些string， swap可以进行优化，这是一个失败的东西：而不是复制，它可以交换所有字符。 我没有尝试过这个，但是我确信这会带来很大的变化。 而且，再说一遍，如果有人不听，基数的大小对运行时间几乎没有影响 – 这意味着我应该尝试实现FryGuy，Stephan和EvilTeach提出的build议。
嗯，是的，顺便说一下： caching局部性是一个值得注意的因素 ：从1Mstring开始，运行时不再线性增加。 但是，这可以很容易地修复：我使用插入sorting为小的子集（<= 20个string） – 而不是随机黑客build议的mergesort。 显然，这样的performance甚至比mergesort这样的小列表更好（见我链接的第一篇论文）。
namespace seqan { template <typename It, typename F, typename T> inline void prescan(It front, It back, F op, T const& id) { using namespace std; if (front == back) return; typename iterator_traits<It>::value_type accu = *front; *front++ = id; for (; front != back; ++front) { swap(*front, accu); accu = op(accu, *front); } } template <typename TIter, typename TSize, unsigned int RADIX> inline void radix_permute(TIter front, TIter back, TSize (& bounds)[RADIX], TSize base) { for (TIter i = front; i != back; ++i) ++bounds[static_cast<unsigned int>((*i)[base])]; TSize fronts[RADIX]; std::copy(bounds, bounds + RADIX, fronts); prescan(fronts, fronts + RADIX, std::plus<TSize>(), 0); std::transform(bounds, bounds + RADIX, fronts, bounds, plus<TSize>()); TSize active_base = 0; for (TIter i = front; i != back; ) { if (active_base == RADIX - 1) return; while (fronts[active_base] >= bounds[active_base]) if (++active_base == RADIX - 1) return; TSize current_base = static_cast<unsigned int>((*i)[base]); if (current_base <= active_base) ++i; else std::iter_swap(i, front + fronts[current_base]); ++fronts[current_base]; } } template <typename TIter, typename TSize> inline void insertion_sort(TIter front, TIter back, TSize base) { typedef typename Value<TIter>::Type T; struct { TSize base, len; bool operator ()(T const& a, T const& b) { for (TSize i = base; i < len; ++i) if (a[i] < b[i]) return true; else if (a[i] > b[i]) return false; return false; } } cmp = { base, length(*front) }; // No closures yet. :-( for (TIter i = front + 1; i != back; ++i) { T value = *i; TIter j = i; for ( ; j != front && cmp(value, *(j - 1)); --j) *j = *(j - 1); if (j != i) *j = value; } } template <typename TIter, typename TSize, unsigned int RADIX> inline void radix(TIter top, TIter front, TIter back, TSize base, TSize (& parent_bounds)[RADIX], TSize next) { if (back - front > 20) { TSize bounds[RADIX] = { 0 }; radix_permute(front, back, bounds, base); // Sort current bucket recursively by suffix. if (base < length(*front) - 1) radix(front, front, front + bounds[0], base + 1, bounds, static_cast<TSize>(0)); } else if (back - front > 1) insertion_sort(front, back, base); // Sort next buckets on same level recursively. if (next == RADIX - 1) return; radix(top, top + parent_bounds[next], top + parent_bounds[next + 1], base, parent_bounds, next + 1); } template <typename TIter> inline void radix_sort(TIter front, TIter back) { typedef typename Container<TIter>::Type TStringSet; typedef typename Value<TStringSet>::Type TString; typedef typename Value<TString>::Type TChar; typedef typename Size<TStringSet>::Type TSize; TSize const RADIX = ValueSize<TChar>::VALUE; TSize bounds[RADIX]; radix(front, front, back, static_cast<TSize>(0), bounds, RADIX - 1); } } // namespace seqan
你可以通过以位为单位编码序列来降低内存需求。 您正在查看排列，因此，对于长度为2的“ACGT”，即16个州或4位。 对于长度为3的64个状态，可以用6位编码。 因此，对于序列中的每个字母，看起来像是2位，或者像你说的那样，对于16个字符，大约是32位。
如果有办法减less有效的“单词”的数量，进一步的压缩可能是可能的。
因此，对于长度为3的序列，可以创build64个桶，可能大小为uint32或uint64。 将它们初始化为零。 迭代3个字符序列的非常大的列表，并按照上面的方式对它们进行编码。 使用这个作为下标，并增加该桶。
重复这个操作，直到所有的序列都被处理完毕。
接下来，重新生成你的列表。
为了在该桶中find的计数顺序遍历64个桶，生成由该桶表示的序列的多个实例。
当所有的桶都被迭代后，你有你的sorting数组。
一个4的序列，增加2位，所以会有256个桶。 一个5的序列，增加2位，所以会有1024个桶。
在某个时候，水桶的数量将接近你的极限。 如果您从文件中读取序列，而不是将其保存在内存中，则可以使用更多内存用于存储桶。
我认为这会比原地进行更快，因为水桶很可能适合你的工作。
这是一个显示技术的黑客
#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
我也会试验分组到更多的桶中，例如，如果你的string是：
GATTACA
第一个MSB调用将返回GATT桶（256个总桶），这样你减less了基于磁盘的缓冲区的分支。 这可能会也可能不会提高性能，所以试试吧。
我要出去，build议你切换到堆/ heapsort实现。 这个build议带有一些假设：
- 您可以控制数据的读取
- 只要“开始”sorting，就可以对sorting后的数据做一些有意义的事情。
堆/堆sorting的美妙之处在于，您可以在读取数据时构build堆，并且可以在构build堆的那一刻开始获取结果。
让我们退后一步。 如果你非常幸运，你可以asynchronous地读取数据（也就是说，你可以发布某种读取请求，并在一些数据准备就绪时收到通知），然后你可以build立一堆堆，而你正在等待下一个数据块 – 即使从磁盘。 通常情况下，这种方法可以在获取数据所花费的时间之后，将大部分的sorting成本都埋在这里。
读取数据后，第一个元素已经可用。 根据您发送数据的位置，这可能会很好。 如果您将其发送给另一个asynchronous阅读器，或者某个并行“事件”模型或UI，则可以随时发送块和块。
这就是说 – 如果你无法控制数据的读取方式，并且它是同步读取的，而且在sorting后的数据完全写出之前没有用处 – 忽略所有这些。 🙁
请参阅维基百科文章：
- 堆sorting
- 二进制堆
性能明智的，你可能想看看更一般的string比较sortingalgorithm。
目前你接触每一个string的每一个元素，但你可以做的更好！
特别是突发sorting非常适合这种情况。 作为一个奖励，因为脉冲串是基于尝试的，所以对于在DNA / RNA中使用的小字母大小来说，它的工作方式非常好，因为你不需要build立任何三元search节点，哈希或其他节点压缩scheme实施。 这些尝试可能对你的后缀数组最终目标也是有用的。
一个体面的通用目的的实现可以在http://sourceforge.net/projects/burstsort/上的source forge上find – 但是它并不在位。
为了比较的目的， http: //www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf中的C-burstsort实现比一些典型工作负载的快速sorting和基数sorting快4-5倍。
你可以尝试使用一个trie 。 对数据sorting只是遍历数据集并插入; 结构是自然sorting的，你可以把它看作类似于B-Tree（除了进行比较，你总是使用指针间接指针）。
caching行为将有利于所有的内部节点，所以你可能不会改善; 但是你也可以摆弄你的trie的分支因子（确保每个节点都适合一个高速caching行，将类似堆的trie节点分配为表示水平顺序遍历的连续数组）。 由于尝试也是数字结构（对于长度为k的元素，O（k）插入/查找/删除），您应该具有基数sorting的有竞争力的性能。
我会分解string的压缩比特表示。 Burstsort声称比基数types有更好的地方性，保持额外的空间使用率，以突破尝试代替古典尝试。 原始纸张有测量。
您将需要看看Drs 的大规模基因组序列处理 。 笠原和森下。
由四个核苷酸字母A，C，G和T组成的string可以被专门编码为整数，以便更快地处理。 基数sorting是书中讨论的许多algorithm之一; 你应该能够适应这个问题的接受的答案，并看到一个很大的性能改善。
“ 没有额外的空间基数sorting ”是一个解决您的问题的论文。
基数sorting不是caching意识，不是大集合最快的sortingalgorithm。 你可以看看：
- ti7qsort 。 ti7qsort是最快的整数sorting（可以用于小尺寸的string）。
- 内联QSORT
- stringsorting
您也可以使用压缩和编码您的DNA的每个字母2位之前存储到sorting数组。
dsimcha的MSB基数sorting看起来不错，但是Nils更接近问题的核心，观察到caching局部性是在大问题大小的情况下会造成什么结果。
我build议一个非常简单的方法：
- 根据经验估计基数sorting的最大大小
m。 - 一次读取
m元素的块，按照基数对它们进行sorting，然后将它们写出（如果内存足够，则将其写入内存缓冲区，否则将被写入）直到耗尽input。 - 合并sorting的结果块。
Mergesort是我知道的最容易caching的sortingalgorithm：“读取数组A或B的下一个项目，然后写入一个项目到输出缓冲区。 它在磁带驱动器上高效运行。 它确实需要2n空间来分类n项目，但是我敢打赌，你将会看到大大改进的caching局部性将会变得不那么重要 – 如果你使用的是非原位基数sorting，那么你需要额外的空间无论如何。
最后请注意，mergesort可以在不recursion的情况下实现，实际上这样做可以清楚地表明真正的线性内存访问模式。
看起来你已经解决了这个问题，但是为了logging，似乎可用的就地基数sorting的一个版本是“美国国旗sorting”。 这里描述： 工程基数sorting 。 总体思路是对每个angular色进行2次传递 – 首先统计每个angular色的数量，这样就可以将input数组细分为分组。 然后再次通过，交换每个元素到正确的箱子。 现在recursion地sorting下一个字符位置上的每个bin。
首先，考虑你的问题的编码。 摆脱string，用二进制表示replace它们。 使用第一个字节来表示长度+编码。 或者，在四字节边界处使用固定长度表示法。 然后基数sorting变得更容易。 对于基数sorting来说，最重要的是在内循环的热点处没有exception处理。
好的，我想了解一下四元问题。 你需要像朱迪树这样的解决scheme。 下一个解决scheme可以处理变长string; 对于固定的长度只是删除长度位，这实际上更容易。
分配16个指针的块。 指针的最低有效位可以被重用，因为你的块总是alignment的。 您可能需要一个特殊的存储分配器（将大容量存储分成更小的块）。 有几种不同types的块：
- 使用7位长度可变的string进行编码。 当他们填满时，你把它们replace成：
- 位置对接下来的两个字符进行编码，您有16个指向下一个块的指针，结束于：
- 位图编码string的最后三个字符。
对于每种块，您需要在LSB中存储不同的信息。 由于你有可变长度的string，所以你也需要存储string结尾，最后一种块只能用于最长的string。 当你深入到结构中时，7个长度的位应该被replace为更less。
这为您提供了一个合理的快速和非常有效的内存有序的string存储。 它会performance得有点像一个线索 。 要做到这一点，请确保构build足够的unit testing。 你想覆盖所有的块转换。 你只想从第二种块开始。
为了获得更高的性能，您可能需要添加不同的块types和更大的块。 如果块总是相同的大小和足够大，你可以使用更less的指针位。 对于16个指针的块大小，在32位地址空间中已经有一个空闲字节。 看看有趣的块types的朱迪树文档。 基本上，您需要为空间（和运行时间）权衡添加代码和工程时间
您可能希望从头四个字符的256宽直接基数开始。 这提供了一个体面的空间/时间折衷。 In this implementation, you get much less memory overhead than with a simple trie; it is approximately three times smaller (I haven't measured). O(n) is no problem if the constant is low enough, as you noticed when comparing with the O(n log n) quicksort.
Are you interested in handling doubles? With short sequences, there are going to be. Adapting the blocks to handle counts is tricky, but it can be very space-efficient.
