205

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

予備

文字「A」、「C」、「G」、「T」のみを使用する小さな固定長の文字列が大量にあり(そう、ご想像のとおり: DNA )、それらを並べ替えます。

現時点では、 STLのすべての一般的な実装で、std::sortwhich uses introsortを使用しています。これは非常にうまく機能します。しかし、基数ソートは問題セットに完全に適合し、実際にははるかにうまく機能するはずだと確信しています。

詳細

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

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

ユースケース

理想的には、このアルゴリズムは、DNA と DNA5 (追加のワイルドカード文字「N」を許可する)、またはIUPAC 曖昧コードを含む DNA (結果として 16 の異なる値) に対して、2 から 100 の間の任意の文字列長で動作する必要があります。ただし、これらすべてのケースをカバーすることはできないことを認識しているため、速度が向上したことに満足しています。コードは、どのアルゴリズムにディスパッチするかを動的に決定できます。

リサーチ

残念ながら、基数ソートに関するウィキペディアの記事は役に立ちません。インプレース バリアントに関するセクションは完全にゴミです。基数ソートに関するNIST-DADS セクションはほとんど存在しません。アルゴリズム「MSL」を説明するEfficient Adaptive In-Place Radix Sortingという有望な論文があります。残念ながら、この論文も残念です。

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

まず、アルゴリズムにはいくつかの誤りがあり、多くのことが説明されていません。特に、再帰呼び出しについては詳しく説明していません (現在のシフト値とマスク値を計算するために、ポインターをインクリメントまたは削減すると仮定しています)。また、関数dest_groupanddest_addressを定義せずに使用します。これらを効率的に実装する方法がわかりません (つまり、O(1) では、少なくともdest_address自明ではありません)。

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

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

要約すると、動作する参照実装、または少なくとも DNA 文字列で動作する動作中の基数ソートの適切な疑似コード/説明を見つける希望はありますか?

4

15 に答える 15

62

さて、これは 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 倍高速です。

于 2009-01-20T21:19:25.870 に答える
21

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

理由：

並べ替えは入力配列に対して線形読み取りを行いますが、すべての書き込みはほぼランダムになります。特定のNから上に向かって、これは書き込みごとのキャッシュミスに要約されます。このキャッシュミスは、アルゴリズムの速度を低下させる原因です。所定の位置にあるかどうかにかかわらず、この効果は変わりません。

これでは質問に直接答えられないことはわかっていますが、並べ替えがボトルネックである場合は、前処理ステップとして並べ替えアルゴリズムの近くを確認することをお勧めします（ソフトヒープのwikiページで開始できる場合があります）。

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

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

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

于 2009-01-20T21:36:30.173 に答える
9

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

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

したがって、長さ 3 のシーケンスの場合、uint32 または uint64 のサイズの 64 個のバケットを作成できます。それらをゼロに初期化します。3 つの char シーケンスの非常に大きなリストを繰り返し処理し、上記のようにエンコードします。これを添え字として使用し、そのバケットを増やします。
すべてのシーケンスが処理されるまで、これを繰り返します。

次に、リストを再生成します。

64 個のバケットを順番に繰り返し処理し、そのバケットで見つかったカウントについて、そのバケットによって表されるシーケンスのインスタンスをその数だけ生成します。
すべてのバケットが繰り返されると、並べ替えられた配列が得られます。

4 のシーケンスでは 2 ビットが追加されるため、256 バケットになります。5 のシーケンスは 2 ビットを追加するため、1024 バケットになります。

ある時点で、バケットの数が制限に近づきます。シーケンスをメモリに保持するのではなく、ファイルから読み取ると、バケットに使用できるメモリが増えます。

バケットは作業セット内に収まる可能性が高いため、in situ で並べ替えを行うよりも高速になると思います。

これはテクニックを示すハックです

#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;
}
于 2009-01-20T22:48:02.010 に答える
6

データセットが非常に大きい場合は、ディスクベースのバッファアプローチが最適だと思います。

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バケット）を返します。これにより、ディスクベースのバッファーのブランチを少なくすることができます。これによりパフォーマンスが向上する場合と向上しない場合があるため、試してみてください。

于 2009-01-20T21:24:34.970 に答える
6

手足に出て、ヒープ/ヒープソートの実装に切り替えることをお勧めします。この提案には、いくつかの前提があります。

  1. データの読み取りを制御します
  2. ソートされたデータを「開始」するとすぐに、ソートされたデータで何か意味のあることを行うことができます。

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

一歩下がってみましょう。幸運にもデータを非同期で読み取ることができる場合（つまり、ある種の読み取り要求を投稿して、データの準備ができたときに通知を受け取ることができる場合）、待機中にヒープのチャンクを構築できます。入ってくるデータの次のチャンク-ディスクからでも。多くの場合、このアプローチは、データの取得に費やされた時間の背後に、ソートの半分のコストの大部分を埋めることができます。

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

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

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

于 2009-01-20T22:00:44.570 に答える
5

「余分なスペースのない基数の並べ替え」は、あなたの問題に対処する論文です。

于 2010-08-05T07:07:06.063 に答える
4

DrsによるLarge-scale Genome Sequence Processingをご覧ください。笠原と森下。

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

于 2010-01-23T18:17:44.187 に答える
3

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

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

于 2009-01-21T05:05:16.630 に答える
3

文字列のパックビット表現をバーストソートします。Burstsort は、基数ソートよりも局所性がはるかに優れていると主張されており、従来の試行の代わりにバースト試行を使用することで、余分なスペースの使用を抑えることができます。元の紙には寸法があります。

于 2009-01-24T22:11:30.693 に答える
2

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

于 2009-01-23T23:50:35.627 に答える
2

基数ソートはキャッシュを意識しておらず、大規模なセットの最速のソート アルゴリズムではありません。あなたは見ることができます：

また、並べ替え配列に格納する前に、圧縮を使用して DNA の各文字を 2 ビットにエンコードすることもできます。

于 2009-06-14T10:25:39.277 に答える
1

dsimcha の MSB 基数ソートは良さそうに見えますが、Nils は問題の核心に迫っており、問題のサイズが大きくなるとキャッシュの局所性が問題になるという観察結果が出ています。

私は非常に単純なアプローチを提案します：

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

Mergesort は、私が知っている中で最もキャッシュに適した並べ替えアルゴリズムです。「配列 A または B から次の項目を読み取り、出力バッファーに項目を書き込みます。」テープ ドライブ上で効率的に動作します。2nアイテムを並べ替えるにはスペースが必要ですnが、キャッシュの局所性が大幅に改善されたことで、それが重要でなくなることは間違いありません。

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

于 2009-01-21T11:40:03.013 に答える
1

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

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

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

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

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

これにより、ソートされた文字列のかなり高速でメモリ効率の高いストレージが提供されます。これは、 triのように振る舞います。これを機能させるには、十分な単体テストを作成してください。すべてのブロック遷移をカバーしたい。2 番目の種類のブロックのみから開始します。

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

最初の 4 文字を 256 幅の直接基数で開始することをお勧めします。これにより、適切な空間/時間のトレードオフが提供されます。この実装では、単純な試行よりもはるかに少ないメモリ オーバーヘッドが得られます。約 3 分の 1 です (測定していません)。O(n log n) クイックソートと比較したときに気づいたように、定数が十分に低い場合、O(n) は問題ありません。

ダブルスの扱いに興味がありますか？短いシーケンスでは、そうなるでしょう。カウントを処理するようにブロックを適応させるのは難しいですが、非常にスペース効率が良い場合があります。

于 2009-01-20T22:45:21.907 に答える