webentwicklung-frage-antwort-db.com.de

In-Place-Radix-Sortierung

Dies ist ein langer Text. Bitte bei mir tragen. Auf den Punkt gebracht lautet die Frage: Gibt es einen funktionsfähigen Radix-Sortieralgorithmus ?


Vorläufig

Ich habe eine große Anzahl von kleinen Strings mit fester Länge, die nur die Buchstaben "A", "C", "G" und "T" verwenden (ja, Sie haben es erraten) : DNA ), die ich sortieren möchte.

Im Moment benutze ich std::sort mit introsort in allen gängigen Implementierungen des STL . Das funktioniert ganz gut. Ich bin jedoch überzeugt, dass radix sort perfekt zu meinem Problem passt und in der Praxis viel besser funktionieren sollte.

Einzelheiten

Ich habe diese Annahme mit einer sehr naiven Implementierung getestet, und für relativ kleine Eingaben (in der Größenordnung von 10.000) stimmte dies (zumindest mehr als doppelt so schnell). Die Laufzeit nimmt jedoch ab, wenn das Problem größer wird ( [~ # ~] n [~ # ~]> 5.000.000).

Der Grund liegt auf der Hand: radix sort erfordert das Kopieren der gesamten Daten (tatsächlich mehr als einmal in meiner naiven Implementierung). Dies bedeutet, dass ich ~ 4 GiB in meinen Hauptspeicher gesteckt habe, was offensichtlich die Leistung beeinträchtigt. Selbst wenn dies nicht der Fall wäre, kann ich es mir nicht leisten, so viel Speicher zu verwenden, da das Problem tatsächlich so groß ist noch größer werden.

Anwendungsfälle

Im Idealfall sollte dieser Algorithmus mit einer beliebigen Stringlänge zwischen 2 und 100 funktionieren, sowohl für DNA als auch für DNA5 (die ein zusätzliches Platzhalterzeichen "N" zulässt) oder sogar für DNA mit IUPAC Mehrdeutigkeitscodes (was zu 16 unterschiedlichen Werten führt). Mir ist jedoch klar, dass all diese Fälle nicht abgedeckt werden können, und ich bin mit jeder Geschwindigkeitsverbesserung zufrieden, die ich erhalte. Der Code kann dynamisch entscheiden, an welchen Algorithmus gesendet werden soll.

Forschung

Leider ist der Wikipedia-Artikel über radix sort unbrauchbar. Der Abschnitt über eine In-Place-Variante ist kompletter Müll. Der NIST-DADS-Abschnitt über die Radix-Sortierung ist so gut wie nicht vorhanden. Es gibt ein vielversprechend klingendes Papier namens Efficient Adaptive In-Place Radix Sorting , das den Algorithmus „MSL“ beschreibt. Leider ist auch dieses Papier enttäuschend.

Insbesondere gibt es die folgenden Dinge.

Erstens enthält der Algorithmus mehrere Fehler und lässt vieles unerklärt. Insbesondere wird der Rekursionsaufruf nicht detailliert dargestellt (ich gehe einfach davon aus, dass er einen Zeiger erhöht oder verringert, um die aktuellen Verschiebungs- und Maskenwerte zu berechnen). Außerdem werden die Funktionen dest_group und dest_address ohne Angabe von Definitionen. Ich verstehe nicht, wie ich diese effizient umsetzen kann (dh in O (1); zumindest dest_address ist nicht trivial).

Last but not least erreicht der Algorithmus In-Place-Ness, indem Array-Indizes mit Elementen innerhalb des Input-Arrays ausgetauscht werden. Dies funktioniert offensichtlich nur bei numerischen Arrays. Ich muss es für Saiten verwenden. Natürlich könnte ich auch einfach mit starker Tipparbeit weitermachen und davon ausgehen, dass der Speicher es toleriert, einen Index dort zu speichern, wo er nicht hingehört. Dies funktioniert jedoch nur, solange ich meine Zeichenfolgen in 32-Bit-Speicher komprimieren kann (32-Bit-Ganzzahlen vorausgesetzt). Das sind nur 16 Zeichen (ignorieren wir für den Moment, dass 16> log (5.000.000)).

Ein anderes Papier von einem der Autoren gibt keine genaue Beschreibung, aber es gibt die Laufzeit von MSL als sublinear an, was absolut falsch ist.

Um zusammenzufassen: Besteht die Hoffnung, eine funktionierende Referenzimplementierung oder zumindest einen guten Pseudocode/eine Beschreibung einer funktionsfähigen Radix-Sortierung zu finden, die mit DNA-Strings funktioniert?

193
Konrad Rudolph

Nun, hier ist eine einfache Implementierung einer MSD-Radix-Sortierung für DNA. Es ist in D geschrieben, da dies die Sprache ist, die ich am häufigsten verwende und daher am seltensten alberne Fehler mache, aber es könnte leicht in eine andere Sprache übersetzt werden. Es ist vorhanden, erfordert aber 2 * seq.length durchläuft das Array.

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);
   }
}

Offensichtlich ist dies spezifisch für DNA, im Gegensatz zu allgemein, aber es sollte schnell sein.

Bearbeiten:

Ich wurde neugierig, ob dieser Code tatsächlich funktioniert, also habe ich ihn getestet/debuggt und darauf gewartet, dass mein eigener Bioinformatik-Code ausgeführt wird. Die obige Version ist jetzt tatsächlich getestet und funktioniert. Bei 10 Millionen Sequenzen mit jeweils 5 Basen ist dies etwa dreimal schneller als bei einem optimierten Introsort.

58
dsimcha

Ich habe noch nie eine in-place-Radix-Sortierung gesehen, und aufgrund der Art der Radix-Sortierung bezweifle ich, dass sie viel schneller ist als eine in-place-Sortierung, solange das temporäre Array in den Speicher passt.

Grund:

Bei der Sortierung wird das Eingabearray linear gelesen, aber alle Schreibvorgänge erfolgen nahezu zufällig. Ab einem bestimmten N führt dies zu einem Cache-Fehler pro Schreibvorgang. Dieser Cache-Fehler verlangsamt Ihren Algorithmus. Ob es vorhanden ist oder nicht, wird diesen Effekt nicht ändern.

Ich weiß, dass dies Ihre Frage nicht direkt beantworten wird, aber wenn das Sortieren ein Engpass ist, sollten Sie sich Near-Sorting Algorithmen als Vorverarbeitungsschritt (Wiki- Seite auf dem Soft-Heap kann Ihnen den Einstieg erleichtern).

Das könnte einen sehr netten Cache-Lokalitätsschub geben. Eine ausserbörsliche Radix-Sortierung für Lehrbücher erzielt dann eine bessere Leistung. Die Schreibvorgänge werden immer noch nahezu zufällig sein, aber zumindest gruppieren sie sich um dieselben Speicherblöcke und erhöhen so die Cache-Trefferquote.

Ich habe jedoch keine Ahnung, ob es in der Praxis funktioniert.

Übrigens: Wenn Sie sich nur mit DNA-Strings beschäftigen: Sie können ein Zeichen in zwei Bits komprimieren und Ihre Daten ziemlich oft packen. Dies verringert den Speicherbedarf um den Faktor vier gegenüber einer Naiive-Darstellung. Die Adressierung wird komplexer, aber die ALU Ihrer CPU hat ohnehin viel Zeit für alle Cache-Misses.

20

Sie können den Speicherbedarf mit Sicherheit verringern, indem Sie die Sequenz in Bits codieren. Sie betrachten Permutationen also für die Länge 2 mit "ACGT", das sind 16 Zustände oder 4 Bits. Für Länge 3 sind das 64 Zustände, die in 6 Bits codiert werden können. So sieht es aus wie 2 Bits für jeden Buchstaben in der Sequenz oder ungefähr 32 Bits für 16 Zeichen, wie Sie sagten.

Wenn es eine Möglichkeit gibt, die Anzahl der gültigen "Wörter" zu verringern, ist möglicherweise eine weitere Komprimierung möglich.

Für Sequenzen der Länge 3 könnten also 64 Buckets erstellt werden, möglicherweise mit der Größe uint32 oder uint64. Initialisieren Sie sie auf Null. Durchlaufen Sie Ihre sehr große Liste von 3 Zeichensequenzen und codieren Sie sie wie oben beschrieben. Verwenden Sie dies als Index und erhöhen Sie diesen Bucket.
Wiederholen Sie diesen Vorgang, bis alle Ihre Sequenzen verarbeitet wurden.

Generieren Sie als Nächstes Ihre Liste neu.

Durchlaufen Sie die 64 Buckets, um für die in diesem Bucket gefundene Anzahl so viele Instanzen der durch diesen Bucket dargestellten Sequenz zu generieren.
Wenn alle Buckets iteriert wurden, haben Sie Ihr sortiertes Array.

Bei einer Folge von 4 werden 2 Bits hinzugefügt, sodass 256 Buckets vorhanden sind. Bei einer Folge von 5 werden 2 Bits hinzugefügt, sodass 1024 Buckets vorhanden sind.

Irgendwann nähert sich die Anzahl der Eimer Ihren Grenzen. Wenn Sie die Sequenzen aus einer Datei lesen, anstatt sie im Speicher zu belassen, steht mehr Speicher für Buckets zur Verfügung.

Ich denke, dies wäre schneller als das Sortieren vor Ort, da die Eimer wahrscheinlich in Ihr Arbeitsset passen.

Hier ist ein Hack, der die Technik zeigt

#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;
}
8
EvilTeach

Wenn Ihre Datenmenge so groß ist, würde ich denken, dass ein datenträgerbasierter Pufferansatz am besten ist:

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

Ich würde auch experimentieren, in eine größere Anzahl von Eimern zu gruppieren, zum Beispiel, wenn Ihre Zeichenfolge war:

GATTACA

beim ersten MSB-Aufruf wird der Bucket für GATT zurückgegeben (insgesamt 256 Buckets). Auf diese Weise werden weniger Verzweigungen des festplattenbasierten Puffers erstellt. Dies kann die Leistung verbessern oder nicht. Probieren Sie es aus.

6
FryGuy

Ich werde auf einen Ast gehen und vorschlagen, dass Sie zu einer Heap/ Heapsort -Implementierung wechseln. Dieser Vorschlag geht von folgenden Annahmen aus:

  1. Sie steuern das Lesen der Daten
  2. Sie können mit den sortierten Daten etwas Sinnvolles tun, sobald Sie anfangen, sie zu sortieren.

Das Schöne an Heap/Heap-Sort ist, dass Sie den Heap erstellen können, während Sie die Daten lesen. Sobald Sie den Heap erstellt haben, können Sie erste Ergebnisse erzielen.

Lass uns zurücktreten. Wenn Sie so viel Glück haben, dass Sie die Daten asynchron lesen können (dh Sie können eine Art Leseanforderung senden und benachrichtigt werden, wenn einige Daten bereit sind), und dann können Sie einen Teil des Heaps aufbauen, während Sie auf das warten Der nächste Datenblock kommt - sogar von der Festplatte. Dieser Ansatz kann häufig den größten Teil der Kosten für die Hälfte Ihrer Sortierung hinter dem Zeitaufwand für das Abrufen der Daten begraben.

Sobald Sie die Daten gelesen haben, ist das erste Element bereits verfügbar. Je nachdem, wohin Sie die Daten senden, kann dies sehr hilfreich sein. Wenn Sie es an einen anderen asynchronen Reader oder an ein paralleles 'Event'-Modell oder eine parallele Benutzeroberfläche senden, können Sie währenddessen Blöcke und Blöcke senden.

Das heißt, wenn Sie keine Kontrolle darüber haben, wie die Daten gelesen werden, und sie synchron gelesen werden und Sie keine Verwendung für die sortierten Daten haben, bis sie vollständig ausgeschrieben sind, ignorieren Sie dies alles. :(

Siehe die Wikipedia-Artikel:

6
Joe

In Bezug auf die Leistung möchten Sie möglicherweise einen allgemeineren Algorithmus zum Sortieren von Zeichenfolgen im Vergleich betrachten.

Momentan berühren Sie jedes Element jeder Saite, aber Sie können es besser machen!

Insbesondere passt ein Burst-Sortierung sehr gut in diesen Fall. Als Bonus, da Burstsort auf Versuchen basiert, funktioniert es lächerlich gut für die kleinen Alphabetgrößen, die in DNA/RNA verwendet werden, da Sie keine Art von ternärem Suchknoten, Hash oder anderem Trie-Knoten-Komprimierungsschema in das Programm einbauen müssen versuchen Umsetzung. Die Versuche können auch für Ihr Endziel nützlich sein, das einem Suffix-Array ähnelt.

Eine ordentliche Allzweckimplementierung von Burstsort ist auf Source Forge unter http://sourceforge.net/projects/burstsort/ verfügbar, ist jedoch nicht vorhanden.

Zu Vergleichszwecken wird die C-Burstsort-Implementierung unter http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf 4-5x schneller als bei quicksort und radix sortiert nach typischen Workloads.

4
Edward KMETT

Schauen Sie sich Large-scale Genome Sequence Processing von Dr. Kasahara und Morishita.

Zeichenfolgen, die aus den vier Nucleotidbuchstaben A, C, G und T bestehen, können speziell in Ganzzahlen codiert werden, um eine viel schnellere Verarbeitung zu erreichen. Die Radix-Sortierung ist unter vielen Algorithmen, die in dem Buch diskutiert werden. Sie sollten in der Lage sein, die akzeptierte Antwort auf diese Frage anzupassen und eine große Leistungsverbesserung zu sehen.

4
Rudiger

" Radix-Sortierung ohne zusätzlichen Speicherplatz " ist ein Artikel, der sich mit Ihrem Problem befasst.

4
eig

Sie könnten versuchen, ein trie zu verwenden. Das Sortieren der Daten erfolgt einfach durch Iteration und Einfügen des Datensatzes. Die Struktur ist natürlich sortiert, und Sie können sich vorstellen, dass sie einem B-Tree ähnelt (außer dass Sie anstelle von Vergleichen immer Zeiger-Indirektionen verwenden). .

Das Caching-Verhalten wird für alle internen Knoten von Vorteil sein, sodass Sie dies wahrscheinlich nicht verbessern werden. Sie können aber auch mit dem Verzweigungsfaktor Ihres Tries experimentieren (stellen Sie sicher, dass jeder Knoten in eine einzelne Cache-Zeile passt, und weisen Sie die einem Heap ähnlichen Trie-Knoten als zusammenhängendes Array zu, das einen Durchlauf in Ebenenreihenfolge darstellt). Da es sich bei Versuchen auch um digitale Strukturen handelt (O (k) Einfügen/Suchen/Löschen für Elemente der Länge k), sollten Sie eine konkurrenzfähige Leistung zu einer radix-Sortierung haben.

3
Tom

Ich würde burstsort eine gepackte Bit-Darstellung der Zeichenfolgen. Es wird behauptet, dass Burstsort eine viel bessere Lokalität hat als Radix-Sortierungen, wodurch der zusätzliche Platzbedarf durch Burst-Versuche anstelle von klassischen Versuchen verringert wird. Das Originalpapier hat Maße.

3
Darius Bacon

Radix-Sort ist nicht cachebewusst und nicht der schnellste Sortieralgorithmus für große Mengen. Sie können sich ansehen:

Sie können auch die Komprimierung verwenden und jeden Buchstaben Ihrer DNA in 2 Bits codieren, bevor Sie sie in das Sortierarray speichern.

2
bill

die MSB-Radix-Sortierung von dsimcha sieht gut aus, aber Nils nähert sich dem Kern des Problems mit der Beobachtung, dass die Cache-Lokalität bei großen Problemgrößen tödlich ist.

Ich schlage einen sehr einfachen Ansatz vor:

  1. Schätzen Sie empirisch die größte Größe m, für die eine Radix-Sortierung effizient ist.
  2. Lesen Sie Blöcke von m Elementen gleichzeitig, sortieren Sie sie nach dem Radix und schreiben Sie sie aus (in einen Speicherpuffer, wenn Sie über genügend Speicher verfügen, ansonsten in eine Datei), bis Sie Ihre Eingabe erschöpfen.
  3. Mergesort die resultierenden sortierten Blöcke.

Mergesort ist der cachefreundlichste Sortieralgorithmus, den ich kenne: "Lies das nächste Element von Array A oder B und schreibe dann ein Element in den Ausgabepuffer." Es läuft effizient auf Bandlaufwerken . Es erfordert 2n Speicherplatz zum Sortieren von n Elementen, aber ich wette, dass die deutlich verbesserte Cache-Lokalität das unwichtig macht - und wenn Sie eine nicht vorhandene Radix-Sortierung verwenden, müssen Sie dies tun dieser zusätzliche Raum sowieso.

Beachten Sie abschließend, dass Mergesort ohne Rekursion implementiert werden kann und auf diese Weise das wahre lineare Speicherzugriffsmuster deutlich wird.

1
j_random_hacker

Denken Sie zunächst an die Kodierung Ihres Problems. Entfernen Sie die Zeichenfolgen und ersetzen Sie sie durch eine binäre Darstellung. Verwenden Sie das erste Byte, um Länge + Codierung anzugeben. Alternativ können Sie eine Darstellung mit fester Länge an einer Vier-Byte-Grenze verwenden. Dann wird die Radix-Sortierung viel einfacher. Für eine Radix-Sortierung ist es am wichtigsten, keine Ausnahmebehandlung am Hotspot der inneren Schleife zu haben.

OK, ich habe ein bisschen mehr über das 4-näre Problem nachgedacht. Sie möchten eine Lösung wie Judy tree dafür. Die nächste Lösung kann Zeichenfolgen mit variabler Länge verarbeiten. Für eine feste Länge entfernen Sie einfach die Längenbits, das macht es tatsächlich einfacher.

Ordnen Sie Blöcke mit 16 Zeigern zu. Das niedrigstwertige Bit der Zeiger kann wiederverwendet werden, da Ihre Blöcke immer ausgerichtet sind. Möglicherweise möchten Sie einen speziellen Speicherzuweiser dafür (Aufteilen eines großen Speichers in kleinere Blöcke). Es gibt verschiedene Arten von Blöcken:

  • Codierung mit 7 Längenbits von Zeichenfolgen variabler Länge. Wenn sie voll sind, ersetzen Sie sie durch:
  • Position kodiert die nächsten zwei Zeichen, Sie haben 16 Zeiger auf die nächsten Blöcke, die mit: enden
  • Bitmap-Codierung der letzten drei Zeichen einer Zeichenfolge.

Für jede Art von Block müssen Sie unterschiedliche Informationen in den LSBs speichern. Da Sie Zeichenfolgen mit variabler Länge haben, müssen Sie auch das Ende der Zeichenfolge speichern, und die letzte Art von Block kann nur für die längsten Zeichenfolgen verwendet werden. Die 7 Längenbits sollten durch weniger ersetzt werden, wenn Sie tiefer in die Struktur vordringen.

Auf diese Weise können sortierte Zeichenfolgen relativ schnell und sehr speichereffizient gespeichert werden. Es wird sich wie ein trie verhalten. Stellen Sie sicher, dass Sie genügend Komponententests erstellen, damit dies funktioniert. Sie möchten alle Blockübergänge erfassen. Sie möchten nur mit der zweiten Art von Block beginnen.

Für noch mehr Leistung möchten Sie möglicherweise verschiedene Blocktypen und eine größere Blockgröße hinzufügen. Wenn die Blöcke immer gleich groß und groß genug sind, können Sie noch weniger Bits für die Zeiger verwenden. Bei einer Blockgröße von 16 Zeigern ist in einem 32-Bit-Adressraum bereits ein Byte frei. In der Judy-Baum-Dokumentation finden Sie interessante Blocktypen. Grundsätzlich fügen Sie Code und Entwicklungszeit für einen Kompromiss zwischen Speicherplatz (und Laufzeit) hinzu

Sie möchten wahrscheinlich mit einem 256-fachen direkten Radix für die ersten vier Zeichen beginnen. Dies bietet einen angemessenen Raum/Zeit-Kompromiss. Bei dieser Implementierung ist der Arbeitsspeicheraufwand wesentlich geringer als bei einem einfachen Versuch. es ist ungefähr dreimal kleiner (ich habe nicht gemessen). O(n) ist kein Problem, wenn die Konstante niedrig genug ist, wie Sie beim Vergleich mit der Quicksorte O (n log n) bemerkt haben.

Interessieren Sie sich für den Umgang mit Doppel? Mit kurzen Sequenzen wird es welche geben. Das Anpassen der Blöcke an die Anzahl ist schwierig, kann jedoch sehr platzsparend sein.

1

Es sieht so aus, als hätten Sie das Problem gelöst, aber im Grunde scheint es, dass eine Version einer funktionsfähigen Radix-Sortierung "American Flag Sort" ist. Es ist hier beschrieben: Engineering Radix Sort . Die allgemeine Idee ist, zwei Durchgänge für jedes Zeichen durchzuführen. Zählen Sie zunächst, wie viele von jedem vorhanden sind, damit Sie das Eingabearray in Klassen unterteilen können. Gehen Sie dann noch einmal durch und tauschen Sie jedes Element in den richtigen Behälter. Sortieren Sie nun rekursiv jedes Fach an der nächsten Zeichenposition.

1
AShelly