就地基数排序

https://stackoverflow.com/questions/463105

19-08-2019
|

题

这是一篇很长的文字。请多多包涵。归根结底，问题是： 是否有可行的就地基数排序算法?

初步的

我有大量的 小定长 仅使用字母“A”、“C”、“G”和“T”的字符串（是的，你已经猜到了：脱氧核糖核酸）我想要排序。

目前，我使用 std::sort 它使用介绍在所有常见的实现中 STL. 。这非常有效。然而，我坚信基数排序完全符合我的问题集并且应该可以工作很多在实践中更好。

细节

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

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

用例

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

研究

不幸的是，维基百科关于基数排序的文章是没用的。关于就地变体的部分完全是垃圾。这 NIST-DADS 关于基数排序的部分几乎不存在。有一篇听起来很有前途的论文叫做高效的自适应就地基数排序它描述了算法“MSL”。不幸的是，这篇论文也令人失望。

具体来说，有以下几件事。

首先，该算法包含多个错误，并且有很多无法解释的地方。特别是，它没有详细说明递归调用（我只是假设它增加或减少一些指针来计算当前的移位和掩码值）。此外，它还使用以下功能 dest_group 和 dest_address 没有给出定义。我不知道如何有效地实现这些（即，在 O(1) 中；至少 dest_address 并非微不足道）。

最后但并非最不重要的一点是，该算法通过将数组索引与输入数组内的元素交换来实现就地性。这显然只适用于数值数组。我需要在字符串上使用它。当然，我可以简单地搞定强类型，然后假设内存能够容忍我将索引存储在不属于它的地方。但这只有当我可以将字符串压缩到 32 位内存中（假设 32 位整数）时才有效。这只有 16 个字符（我们暂时忽略 16 > log(5,000,000)）。

其中一位作者的另一篇论文根本没有给出准确的描述，但它给出了 MSL 的运行时间为次线性，这完全是错误的。

回顾一下:是否有希望找到一个可行的参考实现，或者至少找到一个适用于 DNA 字符串的就地基数排序的良好伪代码/描述？

解决方案

好了，这里有一个简单实现MSD基数排序的DNA中。这是写在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个碱基千万序列，它是关于3×比优化的内省排序更快。

其他提示

我从未见过就地基数排序，并从基数排序的性质我怀疑它比出来的地方排序只要临时数组配合到存储器快得多。

<强>原因：

分选做输入阵列上读出的线性的，但所有的写操作将是几乎随机的。从某个N向上将此归结为每写高速缓存未命中。这缓存未命中就是减慢你的算法。如果是到位与否都不会改变这个效果。

我知道这不会直接回答你的问题，但如果排序是一个瓶颈，你可能想看看近排序算法为前工序（在软堆维基页面可能让你开始）。

这会给人一种很不错的缓存局部性的推动作用。然后课本外的地方基数排序将有更好的表现。写操作仍将接近随机的，但至少他们将围绕相同的内存块，因此提高缓存命中率集群。

我不知道，如果它在实践中工作了，虽然。

顺便说一句：如果你只处理DNA串：您可以压缩焦炭引入两位和包装你的数据了不少。这将减少由系数四个内存要求在naiive表示。解决变得更加复杂，但是你的CPU的ALU有大量的时间在所有的缓存缺失反正花。

可以肯定是由编码比特序列掉落存储要求。您正在寻找的排列是这样，长度2，用“ACGT”这16个州，或4位。对于长度为3，这是64个状态，它可以在6位进行编码。所以，它看起来像2比特序列中的每个字母，或约32个比特为16名字符等你说。

如果有减少的有效“单词”的数量的方式，进一步的压缩是可能的。

因此，对于长度为3的序列，人们可以创建64个水桶，也许尺寸UINT32，或UINT64。它们初始化为零。通过你非常非常大的3个字符序列的列表迭代，并对其进行编码如上。以此为标，并增加该桶。结果直到所有的序列都已经被处理重复此。

接下来，重新生成列表。

通过64轮叶，以便迭代，在该桶中的计数，产生由该桶表示的序列的许多实例。结果当所有的桶都被重复，你有你的有序排列。

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

我也想尝试分组到桶的数量较多，例如，如果你的字符串是：

GATTACA

第一MSB调用将返回关贸总协定（总共256个桶）桶，这样，你让基于磁盘的缓存较少的分支机构。这可能会或可能不会改善性能，所以用它进行试验。

我要出去肢体和建议你换到一堆/heapsort 执行情况。这项建议带有一些假设：

你控制的读取数据
你可以做一些有意义的排序的数据，尽快开始获得排序。

美丽的堆/堆的排序，可以建立堆当你读到的数据，你可以开始取得结果的时候，你已经建立堆。

让我们退后一步。如果你是如此幸运的，你可以读取的异步数据(也就是说，你可以发布某种类型的读请求和通知时，一些数据是准备)，然后你就可以建立一个大块堆在等待下一块的数据来在甚至从磁盘。通常，这种方法可以埋葬的大部分费用的一半的排序在后面的花费的时间得到的数据。

一旦你有读取的数据，第一要素已经可以使用。这取决于你在哪里发送的数据，这可能是巨大的。如果发送到另一个步阅读器，或一些平行的'事件'的模式，或的用户界面，你可以送块和大块你去。

这说的-如果你有没有控制这些数据是如何阅读，这是阅读同时，你有没有用于排序的数据，直到它完全是写出来-忽略所有这一切。:(

看到了维基百科文章：

Heapsort
二进制堆

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

目前，您最终会接触到每根弦的每个元素，但您可以做得更好！

特别是，一个突发排序非常适合这种情况。作为一个额外的好处，由于burstsort是基于尝试的，所以它对于DNA/RNA中使用的小字母大小来说效果非常好，因为你不需要在结构中构建任何类型的三元搜索节点、散列或其他trie节点压缩方案。特里树的实现。这些尝试也可能对您的后缀数组式最终目标有用。

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

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

您会想看看大型基因组序列处理由Drs。笠原和森下。

由四种核苷酸字母串A，C，G和T可特殊编码成整数为多更快的处理。基数排序是在书中讨论了许多算法中;你应该能够接受的答案适应这个问题，并看到一个大的性能提升。

“麦冬，没有额外的空间排序”是纸寻址您的问题

您可以尝试使用特里。排序数据简单地通过数据集迭代并将其插入;该结构自然排序，并可以认为它是类似于B树（除了代替进行比较，则总是使用指针间接寻址）。

缓存行为，将有利于所有的内部节点，所以你可能不会在改善;但可以与你的字典树的分支因子拨弄以及（保证每个节点适合单个高速缓存行，分配类似于堆字典树节点，作为表示电平序遍历的连续阵列）。由于尝试也是数字的结构（O（k）的插入/查找/删除长度为k的元素），就应该有竞争力的性能，以一个基数排序。

我会 burstsort 填充位串的表示。 Burstsort声称比基数排序更好的地方，连爆尝试保持额外的空间使用率下降到位经典的尝试。原始纸具有测量。

基排序不是高速缓存二氧并不是最快的排序的算法对大型集。你可以看看:

ti7qsort.ti7qsort是最快的排序，对于整数(可以用于小固定的尺寸strings)。
内联时，qsort
串的排序

你也可以使用压缩和编码的每个字母你的DNA成2位之前存入的排序阵列。

dsimcha的MSB基数排序看起来不错，但是尼尔斯到靠近心脏的问题的观点，高速缓冲地方是什么杀了你在大问题的大小。

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

根据经验估计的最大尺寸 m 对于这一基数排序是有效的。
读块 m 元素在一个时，基数对它们进行排序，并把它们写出来(到缓冲存储器，如果你有足够的存储器，但无文件)，直到废你的输入。
成了归并排序 所得到的排序块。

成了归并排序是最高速缓存友好的排序的算法我意识到:"阅读的下一个项目从阵A或B，然后写一项目的产出的缓冲器。" 它运行高效率地上 磁带驱动器.它的确需要 2n 空间进行排序 n 的项目，但是我的赌注是多少-改进高速缓冲地，你会看到的将会使这个不重要的--并且如果使用的非就地基数排序，则需要额外的空间。

请注意，最终成了归并排序，可以不用递归，而事实上这样做，清楚表明真正的直线存取模式。

看起来你已经解决了这个问题，但备案，似乎一个可行就地基数排序的一个版本是“美国国旗排序”。排序工程基数：它在这里描述。总的想法是做2遍每个人物 - 第一算多少每次都有的，所以你可以细分输入数组箱。然后通过再次，交换的每个元素到正确的垃圾桶。现在递归排序对下一个字符的位置每个仓。

第一，认为有关的编码的问题。摆脱弦，代替他们由二进制代表性。使用的首字节以表明长编码。或者，使用固定长度表示在四个字节的边界。然后在基数的排序变得更加容易。对于一个基排序，最重要的事情是没有的例外处理在热点的内部循环。

OK，我认为多一点的关于4-没问题。你想要一个解决方案喜欢朱迪树这一点。下一个解决方案可以处理变长串；对于固定长度只是删除长位，这实际上使它更加容易。

分配区块的16个指针。最重要的位的指针可以重复使用，因为你总是会对准。你可能会想要一个特殊的存储器为它(打破了大储存成小块)。有许多各种不同的模块：

编码与7长位的变长串。因为他们填补了，你替代他们通过：
位置编码下面的两个字，你有16个指针到下一个街区，结束：
位编码的最后三个字符串。

对于每一种块，你需要储存中的不同信息位.你有变长串需要存储终串过，最后一种框只能用于最长串。7长位应改为较为你得到更深入的结构。

这为你提供一个合理的快速和非常有效的存储记忆的排串。它将表现得有点像一个 trie.得到这个工作，确保建立足够的单元的测试。你想要复盖所有模块的过渡。你要开始只有第二种块。

甚至更多的性能，可能需要添加块不同类型和尺寸较大的区块。如果该区块总是同样大小的和足够大，可以使用甚至更少位数为指针。用一块大小的指针16,你已经有一个字节免费在一个32位地址的空间。看看Judy树的文件有趣的框类型。基本上，你加入代码和工程时间的空间(和运行时间)的交易

你可能要开始一个256宽的直接基于第一个四字。这提供了一个体面的空间/时间的权衡。在这个实现，你会得到多少存储器开销于一个简单的线索;这是大约三倍(我没有测定)。O(n)是没有问题的，如果变是足够低，因为你注意到在比较O(n记录n)快速排序.

您是否有兴趣在处理翻倍?与短的序列，还有要。调整块来处理重要的是棘手的，但它可以非常空间的效率。

许可以下： CC-BY-SA 和归因

不隶属于 StackOverflow