インプレース基数ソート

https://stackoverflow.com/questions/463105

19-08-2019
|

質問

これは長い文章です。ご了承ください。要約すると、質問は次のとおりです。 実行可能なインプレース基数ソートアルゴリズムはありますか?

予備

膨大な数を持っています 小さい固定長 「A」、「C」、「G」、「T」の文字のみを使用する文字列 (はい、ご想像の通り: DNA）を並べ替えたいと思います。

現時点では、私が使用しているのは、 std::sort を使用するイントロソートすべての一般的な実装では、 STL. 。これは非常にうまく機能します。しかし、私はこう確信しています基数ソート私の問題セットに完全に適合しており、うまくいくはずです 多くの 実際にはもっと良くなります。

詳細

私はこの仮定を非常に単純な実装でテストしましたが、比較的小さな入力 (10,000 程度) の場合、これは当てはまりました (まあ、少なくとも 2 倍以上高速です)。ただし、問題のサイズが大きくなると、実行時間は大幅に低下します (N > 5,000,000).

理由は明らかです。基数ソートでは、データ全体をコピーする必要があります (実際、私の単純な実装では複数回)。これは、メインメモリに最大 4 GiB を追加したことを意味し、明らかにパフォーマンスが低下します。たとえそうでなかったとしても、問題のサイズは実際にはさらに大きくなるため、これほど多くのメモリを使用する余裕はありません。

使用例

理想的には、このアルゴリズムは、DNA だけでなく DNA5 (追加のワイルドカード文字「N」を許可します)、さらには 2 から 100 までの任意の文字列長でも機能するはずです IUPAC あいまいなコード (結果として 16 個の異なる値が得られます)。ただし、これらすべてのケースをカバーできるわけではないことは理解しているので、速度が向上することに満足しています。コードは、どのアルゴリズムにディスパッチするかを動的に決定できます。

研究

残念ながら、基数ソートに関するウィキペディアの記事役に立たない。インプレースバリアントに関するセクションは完全にゴミです。の基数ソートに関する NIST-DADS セクション存在しないに等しい。という有望そうな論文がある。効率的な適応型インプレース基数ソートこれはアルゴリズム「MSL」について説明しています。残念ながら、この論文も残念です。

具体的には以下のようなものがあります。

まず、アルゴリズムにはいくつかの間違いが含まれており、多くのことが説明されていません。特に、再帰呼び出しについては詳しく説明しません (現在のシフト値とマスク値を計算するために、何らかのポインターが増加または減少すると単純に考えています)。また、関数を使用します dest_group そして dest_address 定義を与えずに。これらを効率的に (つまり、O(1) で) 実装する方法がわかりません。少なくとも dest_address それは簡単なことではありません）。

最後に重要なことですが、このアルゴリズムは、配列インデックスを入力配列内の要素と交換することによって、インプレース性を実現します。これは明らかに数値配列でのみ機能します。文字列で使用する必要があります。もちろん、強い型付けをダメにして、本来あるべきでない場所にインデックスを保存してもメモリが許容してくれると仮定して続行することもできます。ただし、これは文字列を 32 ビットのメモリ (32 ビット整数を想定) に詰め込むことができる場合にのみ機能します。これはわずか 16 文字です (ここでは、16 > log(5,000,000) であることは無視しましょう)。

著者の一人による別の論文には正確な説明がまったくありませんが、MSL のランタイムがサブリニアであるとされており、これは完全に間違っています。

要点をまとめると:動作するリファレンス実装、または少なくとも DNA 文字列で動作するインプレース基数ソートの適切な擬似コード/記述を見つける希望はありますか?

解決

さて、これは DNA の MSD 基数ソートの簡単な実装です。D で書かれているのは、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 塩基のシーケンスが 1,000 万個ある場合、最適化されたイントロソートよりも約 3 倍高速になります。

他のヒント

インプレース基数ソートを見たことがありません。また、基数ソートの性質から、一時配列がメモリーに収まる限り、アウトオブプレース・ソートよりもはるかに高速であるとは思えません。

理由：

ソートでは入力配列に対して線形読み取りが行われますが、書き込みはすべてほぼランダムになります。特定の N より上では、これは結局、書き込みごとのキャッシュミスになります。このキャッシュミスがアルゴリズムの速度を低下させる原因となります。設置されていてもいなくても、この効果は変わりません。

これがあなたの質問に直接答えられないことは承知していますが、並べ替えがボトルネックになっている場合は、以下を参照してください。 仕分け間近 としてのアルゴリズム 前処理ステップ (ソフトヒープ上の wiki ページから始めることができるかもしれません)。

これにより、キャッシュの局所性が非常に向上する可能性があります。教科書のような場違いな基数ソートのパフォーマンスが向上します。書き込みは依然としてほぼランダムですが、少なくともメモリの同じチャンクの周りに集中するため、キャッシュヒット率が向上します。

実際にそれがうまくいくかどうかはわかりませんが。

ところで:DNA 文字列のみを扱う場合:char を 2 ビットに圧縮して、かなり多くのデータを詰め込むことができます。これにより、単純な表現に比べてメモリ要件が 4 分の 1 に削減されます。アドレス指定はより複雑になりますが、いずれにせよ、CPU の ALU はすべてのキャッシュミスの間に多くの時間を費やします。

シーケンスをビット単位でエンコードすることで、メモリ要件を確実に減らすことができます。順列を調べているので、長さ 2 の場合、「ACGT」では 16 状態、つまり 4 ビットになります。長さ 3 の場合、64 状態になり、6 ビットでエンコードできます。したがって、シーケンス内の各文字に対して2ビット、またはあなたが言ったように16文字に対して約32ビットのように見えます。

有効な「単語」の数を減らす方法がある場合は、さらに圧縮できる可能性があります。

したがって、長さ 3 のシーケンスの場合、サイズが uint32 または uint64 のバケットを 64 個作成できます。それらをゼロに初期化します。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 呼び出しは GATT のバケット (合計 256 バケット) を返します。これにより、ディスクベースのバッファーの分岐が少なくなります。これによりパフォーマンスが向上する場合とそうでない場合があるため、試してみてください。

思い切ってヒープに切り替えることを提案します/ヒープソート実装。この提案にはいくつかの前提条件が含まれています。

データの読み取りを制御するのはあなたです
並べ替えを「開始」するとすぐに、並べ替えられたデータで何か有意義な作業を行うことができます。

ヒープ/ヒープソートの利点は、データの読み取り中にヒープを構築でき、ヒープを構築した瞬間から結果の取得を開始できることです。

後戻りしましょう。幸運にもデータを非同期で読み取ることができれば (つまり、ある種の読み取りリクエストを送信し、データの準備ができたら通知を受け取ることができます)、そのデータの待機中にヒープのチャンクを構築できます。次に受信するデータのチャンク (ディスクからであっても)。多くの場合、このアプローチにより、並べ替えの半分のコストのほとんどを、データの取得に費やされる時間の背後に埋めることができます。

データを読み取ったら、最初の要素はすでに使用可能になっています。データの送信先によっては、これが素晴らしい場合があります。別の非同期リーダー、または並列の「イベント」モデル、または UI に送信する場合は、送信するたびにチャンクを送信できます。

とはいえ、データの読み取り方法を制御できず、データが同期的に読み取られ、完全に書き出されるまでソートされたデータを使用できない場合は、これをすべて無視してください。:(

ウィキペディアの記事を参照してください。

パフォーマンスの観点から、より一般的な文字列比較ソートアルゴリズムを検討することをお勧めします。

現時点では、すべての文字列のすべての要素を触ることになってしまいますが、もっとうまくできるはずです。

特に、バーストソートはこのケースに非常によく当てはまります。おまけに、バーストソートはトライに基づいているため、三分検索ノード、ハッシュ、その他のトライノード圧縮スキームを構築する必要がないため、DNA/RNA で使用される小さなアルファベットサイズに対して驚くほどうまく機能します。実装してみる。この試行は、サフィックス配列のような最終目標にも役立つ可能性があります。

バーストソートの適切な汎用実装は、次のソースフォージで入手できます。 http://sourceforge.net/projects/burstsort/ - しかし、それは適切な場所にありません。

比較のために、C バーストソートの実装については、次の URL で説明しています。 http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf 一部の一般的なワークロードでは、ベンチマークはクイックソートや基数ソートよりも 4 ～ 5 倍高速です。

見てみたいと思うでしょう大規模ゲノム配列処理博士による。笠原さんと森下さん。

4 つのヌクレオチド文字 A、C、G、および T で構成される文字列は、特別に整数にエンコードできます。 多くの より高速な処理。基数ソートは、この本で説明されている多くのアルゴリズムの 1 つです。この質問に対する受け入れられた回答を適用すると、パフォーマンスが大幅に向上するはずです。

"余分なスペースを使用しない基数ソート」はあなたの問題に対処する論文です。

を使用してみてもよいでしょう試してみる. 。データの並べ替えは、データセットを反復処理して挿入するだけです。構造は自然にソートされており、B ツリーに似ていると考えることができます (比較を行う代わりに、 いつも ポインタの間接参照を使用します)。

キャッシュ動作はすべての内部ノードに有利になるため、おそらくこれを改善することはできません。ただし、トライの分岐要素をいじることもできます (すべてのノードが 1 つのキャッシュラインに収まるようにし、レベル順序のトラバーサルを表す連続した配列としてトライノードをヒープと同様に割り当てます)。try もデジタル構造 (長さ k の要素に対する O(k) 挿入/検索/削除) であるため、基数ソートに匹敵するパフォーマンスが得られるはずです。

私は...するだろうバーストソート文字列のパックビット表現。バーストソートは基数ソートよりも局所性がはるかに優れており、従来のトライの代わりにバーストトライを使用することで余分なスペースの使用量を抑えていると主張されています。オリジナルの紙には寸法が記載されています。

Radix-Sort はキャッシュを意識しないため、大規模なセットに対する最速のソートアルゴリズムではありません。以下を参照できます。

ティ7クソート. 。ti7qsort は、整数の最も高速なソートです (小さい固定サイズの文字列に使用できます)。
インライン QSORT
文字列のソート

圧縮を使用して、ソート配列に保存する前に DNA の各文字を 2 ビットにエンコードすることもできます。

dsimcha の MSB 基数ソートは優れているように見えますが、Nils は、問題のサイズが大きい場合にはキャッシュの局所性が問題の原因であると観察し、問題の核心に迫っています。

非常にシンプルなアプローチをお勧めします。

最大サイズを経験的に推定する m 基数ソートが効率的です。
のブロックを読み取ります m 要素を一度に 1 つずつ基数ソートし、入力がなくなるまで (十分なメモリがある場合はメモリバッファに、そうでない場合はファイルに) 書き込みます。
マージソート 結果としてソートされたブロック。

Mergesort は、私が知っている中で最もキャッシュに適した並べ替えアルゴリズムです。「アレイAまたはBのいずれかから次のアイテムを読み、出力バッファーにアイテムを書きます。」効率的に実行されます テープドライブ. 。必要です 2n 並べ替えるスペース n しかし、私の推測では、キャッシュの局所性が大幅に改善されたため、そのことは重要ではなくなるでしょう。また、非インプレース基数ソートを使用していた場合は、とにかく余分なスペースが必要でした。

最後に、マージソートは再帰なしで実装できることに注意してください。実際、この方法で実装すると、真の線形メモリアクセスパターンが明確になります。

問題は解決したように見えますが、記録のために書いておくと、実行可能なインプレース基数ソートのバージョンの 1 つが「アメリカ国旗ソート」であるようです。ここで説明されています: エンジニアリング基数ソート. 。一般的な考え方は、各文字に対して 2 つのパスを実行することです。まず、それぞれのパスの数を数えて、入力配列をビンに分割できるようにします。次に、各要素を正しいビンに入れ替えて再度実行します。次に、各ビンを次の文字位置で再帰的に並べ替えます。

まず、問題のコーディングについて考えてください。文字列を削除し、バイナリ表現に置き換えます。最初のバイトを使用して長さとエンコーディングを示します。あるいは、4 バイト境界で固定長表現を使用します。そうすれば、基数ソートがはるかに簡単になります。基数ソートの場合、最も重要なことは、内部ループのホットスポットで例外処理を行わないことです。

さて、4 進問題についてもう少し考えてみました。次のような解決策が必要ですジュディ・ツリーこのために。次のソリューションでは、可変長の文字列を処理できます。固定長の場合は長さのビットを削除するだけで、実際には簡単になります。

16 ポインターのブロックを割り当てます。ブロックは常に位置合わせされるため、ポインターの最下位ビットは再利用できます。そのために特別なストレージアロケータが必要になる場合があります (大きなストレージを小さなブロックに分割する)。ブロックにはさまざまな種類があります。

7 ビット長の可変長文字列を使用したエンコード。それらがいっぱいになると、次のように置き換えます。
位置は次の 2 文字をエンコードします。次のブロックへのポインターが 16 個あり、次で終わります。
文字列の最後の 3 文字のビットマップエンコーディング。

ブロックの種類ごとに、LSB に異なる情報を格納する必要があります。可変長の文字列がある場合は、文字列の終わりも格納する必要があり、最後の種類のブロックは最長の文字列にのみ使用できます。構造を深く理解するにつれて、7 つの長さのビットをより少ないビットに置き換える必要があります。

これにより、ソートされた文字列を適度に高速でメモリ効率の高いストレージに保存できます。それはある程度次のように動作します試してみる. 。これを機能させるには、十分な単体テストを作成してください。すべてのブロック遷移をカバーする必要があります。2 番目の種類のブロックのみから始めたいとします。

さらにパフォーマンスを向上させるには、さまざまなブロックタイプを追加し、より大きなサイズのブロックを追加することをお勧めします。ブロックが常に同じサイズで十分な大きさであれば、ポインターにさらに少ないビットを使用できます。ブロックサイズが 16 ポインターの場合、32 ビットアドレス空間にはすでに 1 バイトの空きがあります。興味深いブロックタイプについては、Judy ツリーのドキュメントを参照してください。基本的に、スペース (およびランタイム) のトレードオフのために、コードとエンジニアリング時間を追加します。

最初の 4 文字については、256 幅の直接基数から始めるとよいでしょう。これにより、スペースと時間の適切なトレードオフが実現されます。この実装では、単純なトライよりもメモリオーバーヘッドが大幅に少なくなります。約3倍小さくなります（測定はしていません）。O(n log n) クイックソートと比較したときに気づいたように、定数が十分に低い場合、O(n) は問題ありません。

ダブルスの扱いに興味はありますか？短いシーケンスでは、それが起こるでしょう。ブロックを調整してカウントを処理するのは難しいですが、スペース効率は非常に高くなります。

ライセンス： CC-BY-SA と帰属

所属していません StackOverflow