这是一篇很长的文章。请原谅我。归结起来,问题是:是否存在可行的就地基数排序算法?
初步
我有大量固定长度的小字符串,只使用字母“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字符串上的就地基数排序?
您当然可以通过对序列进行位编码来降低内存需求。
你看到的是排列,对于长度为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;
}