Frage

Dies ist ein langer Text. Bitte bei mir tragen. Kocht, ist die Frage:? Gibt es eine praktikable in-place Radixsort Algorithmus

Vorläufige

ich eine große Anzahl von habe kleine feste Länge Strings, die nur mit den Buchstaben „A“, „C“, „G“ und „T“ (ja, Sie haben erraten es:. DNA ), die ich sortieren möchten

Im Moment benutze ich std::sort die Introsort in allen gängigen Implementierungen der STL . Das funktioniert ganz gut. Ich bin jedoch davon überzeugt, dass Radixsort passt mein Problem perfekt eingestellt und sollte arbeiten viel in der Praxis besser.

Details

Ich habe diese Annahme mit einer sehr naiven Implementierung getestet und für relativ kleine Eingänge (in der Größenordnung von 10.000) war dies wahr (na ja, zumindest mehr als doppelt so schnell). Jedoch Laufzeit verschlechtert bodenlos, wenn das Problem größer wird ( N > 5.000.000).

Der Grund ist offensichtlich: Radixsort erfordert das Kopieren der gesamten Daten (mehr als einmal in meine naive Implementierung, tatsächlich). Das bedeutet, dass ich setze ~ 4 gebt in meinen Hauptspeicher, der offensichtlich Leistung tötet. Auch wenn es nicht der Fall war, kann ich nicht leisten, so viel Speicher zu verwenden, da das Problem sogar noch größer geworden Größen.

Use Cases

Idealerweise sollte dieser Algorithmus mit jeder String-Länge arbeiten, zwischen 2 und 100, für die DNA als auch DNA5 (die einen zusätzlichen Platzhalter „N“ erlaubt), oder sogar mit DNA IUPAC Mehrdeutigkeit Codes (was zu 16 verschiedenen Werten). Aber ich merke, dass alle diese Fälle können nicht abgedeckt werden, also bin ich glücklich mit jeder Geschwindigkeit Verbesserung, die ich bekommen. Der Code kann entscheiden, dynamisch welcher Algorithmus zu versenden.

Forschung

Leider ist die Wikipedia-Artikel über Radix ist eine Art nutzlos. Der Abschnitt über eine in-Place-Variante ist völliger Quatsch. Die NIST-DADS Abschnitt auf Radixsort ist neben nicht existent. Es gibt ein viel versprechende klingendes Papier namens Effiziente Adaptive In-Place Radix Sortierung des beschreibt den Algorithmus „MSL“. Leider ist dieses Papier ist auch enttäuschend.

Insbesondere gibt es folgende Dinge.

Zunächst enthält der Algorithmus mehrere Fehler und lässt viel ungeklärt. Insbesondere spielt es keine Detail die Rekursion Anruf (ich einfach davon ausgehen, dass es erhöht oder verringert einige Zeiger der aktuellen Verschiebung und Maskenwerte zu berechnen). Außerdem verwendet er die Funktionen dest_group und dest_address ohne Definitionen zu geben. Ich sehe nicht, wie diese effizient umzusetzen (dh, in O (1), zumindest dest_address ist nicht trivial).

Last but not least, erreicht der Algorithmus in-place-ness von Array-Indizes mit Elementen innerhalb des Eingangs-Array austauschen. Dies funktioniert natürlich nur auf numerische Arrays. Ich brauche es auf Zeichenketten zu verwenden. Natürlich könnte ich nur starke Typisierung Schraube und unter der Annahme, gehen Sie vor, dass der Speicher meines toleriert einen Index zu speichern, wo es nicht gehört. Das funktioniert aber nur so lange, wie ich meine Saiten in 32 Bit Speicher quetschen kann (vorausgesetzt, 32-Bit-Integer). Das ist nur 16 Zeichen (lassen Sie uns für den Moment ignorieren, dass 16> log (5.000.000)).

Ein weiteres Papier von einem der Autoren gibt keine genaue Beschreibung überhaupt, aber es gives MSL Laufzeit als sublinear was falsch ist flach aus.

Zur Erinnerung : Gibt es Hoffnung, eine Arbeitsreferenzimplementierung zu finden oder zumindest eine gute Pseudo-Code / Beschreibung eines Arbeits in-place Radixsort, die auf DNA-Strings arbeiten

War es hilfreich?

Lösung

Nun, hier ist eine einfache Implementierung eines MSD Radixsort für DNA. Es ist in D geschrieben, weil das die Sprache ist, die ich am meisten verwenden, und deshalb bin am wenigsten wahrscheinlich dumme Fehler in machen, aber es könnte leicht in eine andere Sprache übersetzt werden. Es ist an Ort und Stelle, sondern erfordert 2 * seq.length durch das Feld passiert. 

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 diese Art von DNA-spezifischen, im Gegensatz zu allgemeinen zu sein, aber es sollte schnell sein.

Edit:

Ich habe gespannt, ob dieser Code tatsächlich funktioniert, so dass ich getestet / debuggt es während für meine eigene Bioinformatik Code warten, um zu laufen. Die Version oben jetzt tatsächlich getestet und funktioniert. Für 10 Millionen Sequenzen von 5 Basen jedem, es ist etwa 3x schneller als eine optimierte Introsort.

Andere Tipps

Ich habe noch nie einen in-Place Radixsort, und von der Art der Radix-Art Ich bezweifle, dass es viel schneller als ein deplatzierte Art, solange die temporäre Array paßt in den Speicher.

Begründung:

Die Sortierung hat eine linear auf dem Eingabefeld gelesen, aber alle Schreibvorgänge werden fast zufällig. Ab einem bestimmten N kocht aufwärts dies pro Schreibvorgang in eine Cache-Miss nach unten. Diese Cache-Miss ist das, was Sie Ihren Algorithmus verlangsamt. Wenn es an seinem Platz ist oder nicht, wird diese Wirkung nicht ändern.

Ich weiß, dass dies Ihre Frage nicht direkt beantworten, aber wenn das Sortieren ein Engpass ist, kann man einen Blick auf haben will in der Nähe von Sortier- Algorithmen als Vorverarbeitungsschritt ( die Wiki-Seite auf dem weichen Haufen erhalten Sie möglicherweise gestartet).

Das könnte einen sehr schönen Cache-Lokalität Schub geben. Ein Text-Buch out-of-Platz Radixsort wird dann besser. Die Schreibvorgänge werden immer noch fast zufällig sein, aber zumindest werden sie um die gleichen Stücke von Speichern und als solche Erhöhung der Cache-Trefferquote Cluster.

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

Btw: Wenn Sie mit DNA-Strings zu tun hat nur: Sie können ein Zeichen in zwei Bits komprimieren und Ihre Daten sehr viel packen. Dadurch wird der Speicherbedarf um den Faktor vier über eine naiive Darstellung abgeholzt. Die Adressierung wird komplexer, aber die ALU Ihrer CPU hat viel Zeit sowieso bei allen Cache-Misses zu verbringen.

Sie können sicher durch die Codierung der Sequenz in Bits die Speicheranforderungen fallen. Sie betrachten Permutationen so, für Länge 2, mit „ACGT“, die 16 Staaten oder 4 Bits ist. Für Länge 3, ist, dass 64 Staaten, die in 6 Bits codiert werden kann. So sieht es aus wie zwei Bits für jeden Buchstaben in der Reihenfolge, oder etwa 32 Bits für 16 Zeichen wie Sie gesagt haben.

Wenn es eine Möglichkeit gibt, die Anzahl der gültigen ‚Worte‘, eine weitere Verdichtung kann möglich sein.

zu reduzieren

Also für Sequenzen der Länge 3, ein 64 Eimer schaffen könnte, so bemessen, vielleicht Uint32 oder uint64. Initialisieren Sie sie auf Null. Iteration durch Ihre sehr sehr große Liste von drei char-Sequenzen und kodieren sie, wie oben. Verwenden Sie diese als Index und den Eimer erhöhen.
Wiederholen Sie diesen Vorgang, bis alle Ihre Sequenzen verarbeitet worden sind.

Als nächstes regeneriert Ihre Liste.

Iteration durch die 64 Eimer um, für die Zählung in diesem Eimer gefunden, erzeugt, dass viele Fälle der von diesem Eimer Sequenz, dargestellt.
wenn alle der Eimer iteriert wurden, haben Sie Ihre sortierten Array.

Eine Folge von 4, fügt 2 Bits, so würde es 256 Eimer sein. Eine Sequenz von 5, fügt 2 Bits, so gäbe es 1024 Eimer sein.

An einem gewissen Punkt die Anzahl der Schaufeln Ihre Grenzen nähern. Wenn Sie die Sequenzen aus einer Datei lesen, anstatt sie im Speicher zu halten, mehr Speicher wären für Eimer zur Verfügung.

Ich denke, das wäre schneller als die Art in situ zu tun als die Eimer sind wahrscheinlich in Ihrem Arbeitssatz 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;
}

Wenn Ihr Datensatz so groß ist, dann würde ich denken, dass ein Disk-basierten Puffer Ansatz am besten wäre: 

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 experimentieren auch in eine größere Anzahl von Eimern Gruppierung, zum Beispiel, wenn die Zeichenfolge war: 

GATTACA

der erste MSB Anruf würde den Eimer für GATT (256 insgesamt Eimer) zurückkehren, so Sie weniger Zweige der Disk-basierten Puffer machen. Dies kann oder kann nicht die Leistung verbessern, so experimentieren damit.

Ich werde auf einem Bein gehen und schlagen Sie vor, um einen Heap-Schalter / Heapsort Umsetzung. Dieser Vorschlag kommt mit einigen Annahmen:

  1. Sie steuern das Lesen der Daten
  2. Sie können so schnell etwas Sinnvolles mit den sortierten Daten tun, wie Sie auf ‚Start‘ bekommen es sortiert.

Die Schönheit des Haufens / heap-Art ist, dass man den Haufen bauen können, während Sie die Daten lesen, und Sie können beginnen, Resultate in dem Moment bekommen Sie den Heap aufgebaut haben.

Lassen Sie uns einen Schritt zurück. Wenn Sie so glücklich sind, dass Sie die Daten asynchron lesen können (das heißt, können Sie irgendeine Art von Leseanforderung posten und benachrichtigt, wenn einige Daten bereit ist), und dann kann man ein Stück des Haufens bauen, während Sie für die warten nächster Teil der Daten kommen - auch von der Festplatte. Oft kann dieser Ansatz begraben den größten Teil der Kosten für die Hälfte der hinter der Zeit mit dem Sortieren verbrachte die Daten zu erhalten.

Wenn Sie die Daten gelesen haben, ist das erste Element bereits vorhanden. Je nachdem, wo Sie die Daten senden, können diese groß sein. Wenn Sie es an einen anderen asynchronen Leser, oder ein parallel ‚Ereignis‘ Modell oder UI senden, können Sie Stücke und Brocken senden, wie Sie gehen.

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

Sehen Sie Wikipedia-Artikel:

Performance-weise könnten Sie auf einem allgemeineren String-Vergleich Sortieralgorithmen suchen.

Derzeit wickeln Sie jedes Element jeder Saite berühren, aber man kann es besser!

Insbesondere ein platzen sort ist eine sehr gute Passform für diesen Fall. Als Bonus, da burstsort auf Versuchen basiert, funktioniert es lächerlich gut für die kleinen Alphabet Größen in DNA / RNA verwendet, da Sie ternären Suchknoten, Hash oder andere Trie-Knoten Komprimierungsschema keine Art von in die bauen müssen Trie-Implementierung. Die Versuche können für Ihr Suffix-Array-artige Endziel als auch nützlich sein.

Eine anständige Allzweck Umsetzung burstsort ist auf Sourceforge unter http://sourceforge.net/projects / burstsort / -. aber es ist nicht an Ort und Stelle

Zu Vergleichszwecken Die C-burstsort Implementierung unter http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf Benchmarks 4-5x schneller als Quicksort und Radixsort für einige typische Workloads.

„ Radix Sortierung ohne zusätzlichen Platz “ ist ein Papier Ihre Adressierung Problem dar.

Sie könnten versuchen, eine trie . Sortieren der Daten einfach durch den Datensatz wird iteriert und Einsetzen; die Struktur ist natürlich sortiert, und Sie können mit einem B-Baum der es so ähnlich denken (mit Ausnahme statt Vergleiche zu machen, Sie immer Zeiger Indirekt verwenden).

Caching Verhalten wird alle internen Knoten begünstigen, so dass Sie wahrscheinlich auf, dass nicht verbessern wird; Sie können jedoch mit dem Verzweigungsfaktor der Trie-Geige als auch (sicherzustellen, dass jeder Knoten paßt in eine einzige Cache-Zeile, Trie-Knoten zuweisen ähnlich einen Haufen, als eine zusammenhängende Array, das eine Ebene Ordnung Traversal darstellt). Da versucht auch digitale Strukturen (O (k) insert / finden / für Elemente der Länge k löschen), sollten Sie wettbewerbsfähige Leistung zu einem Radixsort haben.

Ich würde burstsort eine gepackte-Bit-Darstellung der Saiten. Burstsort wird behauptet, viel besser als Lokalität Radixsort zu haben, die zusätzliche Raumnutzung niedrig zu halten mit Burst versucht anstelle der klassischen Versuche. Das Originalpapier hat Messungen.

Radix-Sort ist nicht bewusst Cache und ist nicht der schnellste Sortieralgorithmus für große Mengen. Sie können aussehen:

Sie können auch Komprimierung und jeden Buchstaben Ihrer DNA in 2 Bits codieren, bevor sie in die Sortier Array zu speichern.

dsimcha die MSB Radixsort sieht nett aus, aber Nils kommt näher an das Herz des Problems mit der Beobachtung, dass Cache-Ort ist das, was man im Allgemeinen Problemgrößen tötet.

Ich schlage vor, einen sehr einfachen Ansatz:

  1. schätzen Empirisch die größte Größe m, für die eine Radixsort effizient ist.
  2. Lesen Blöcke m Elemente zu einer Zeit, sie sortieren Radix, und schreiben Sie sie aus (in einen Speicherpuffer, wenn Sie genügend Speicher haben, aber ansonsten zur Datei), bis Sie Ihre Eingabe erschöpfen.
  3. Mergesort die resultierenden sortierten Blöcke.

Mergesort ist der Cache freundliche Sortieralgorithmus Ich bin mir dessen bewusst: „Lesen Sie das nächste Element aus entweder Array A oder B, dann auf ein Element in den Ausgangspuffer schreiben.“ Es läuft effizient auf Bandlaufwerke . Es tut 2n Platz benötigt n Elemente zu sortieren, aber meine Wette ist, dass die stark verbesserte Cache-Lokalität, dass unwichtig machen sehen werden wird - und wenn Sie einen nicht-in-place Radixsort verwendet wurden, brauchte man, dass mehr Platz trotzdem.

Bitte beachten Sie schließlich, dass mergesort ohne Rekursion implementiert wird, und in der Tat tut es auf diese Weise klar das wahre lineare Speicherzugriffsmuster macht.

Es sieht aus wie Sie das Problem gelöst haben, aber für die Aufzeichnung, scheint es, dass eine Version eines tragfähigen in-place Radixsort die „American Flag Sort“ ist. Es ist hier beschrieben: Ingenieur Radix Sort . Die allgemeine Idee ist auf jeden Charakter 2 Pässe zu tun - zuerst zählen, wie viele von jedem Sie haben, so dass Sie die Eingabe-Array in Bins unterteilt werden können. Gehen Sie dann wieder durch, jedes Element in den richtigen Behälter austauschen. Nun rekursiv jedes Fach auf der nächsten Zeichenposition sortieren.

Zuerst denken über die Codierung des Problems. Lassen Sie sich von den Saiten befreien, ersetzen Sie sie durch eine binäre Darstellung. Verwenden Sie den ersten Bytelänge + Codierung anzuzeigen. Alternativ dazu verwenden, um eine festgelegte Länge Darstellung in einer Vier-Byte-Grenze. Dann wird der Radixsort viel einfacher. Für eine Radixsort, ist das Wichtigste nicht die Ausnahmebehandlung bei dem Hot-Spot der inneren Schleife hat.

OK, dachte ich ein bisschen mehr über das 4-nären Problem. Sie wollen eine Lösung wie ein Judy Baum für diese. Die nächste Lösung kann Zeichenfolge variabler Länge verarbeiten; für feste Länge entfernen Sie einfach die Länge Bits, das macht es tatsächlich einfacher.

Weisen Blöcke von 16 Zeigern. Das niedrigstwertige Bit der Zeiger kann wiederverwendet werden, da Ihre Blöcke immer ausgerichtet werden. Sie könnten eine spezielle Lagerung Allocator für sie (Zerschlagung großer Speicher in kleinere Blöcke) werden soll. Es gibt eine Reihe von verschiedenen Arten von Blöcken:

  • Codierung mit 7 Länge Bits von Zeichenfolge variabler Länge. Als sie füllen Sie ersetzen Sie sie durch:
  • Position codiert die nächsten zwei Zeichen, Sie haben 16 Zeiger auf die nächsten Blöcke, endend mit:
  • Bitmap-Codierung der letzten drei Zeichen eines Strings.

Für jede Art von Block, müssen Sie verschiedene Informationen in der LSBs speichern. Wie Sie Zeichenfolge variabler Länge haben müssen Sie auch am Ende der Zeichenfolge zu speichern, und die letzte Art von Block kann nur für die längsten Strings verwendet werden. Die 7 Länge Bits sollte durch weniger ersetzt werden, wie Sie tiefer in die Struktur erhalten.

Dies bietet Ihnen eine einigermaßen schnell und sehr speichereffiziente Speicherung von sortierten Ketten. Es verhält sich ein wenig wie ein trie . Um diese Arbeit zu erhalten, stellen Sie sicher, dass genug Unit-Tests zu bauen. Sie wollen Abdeckung aller Satzübergänge. Sie wollen nur mit der zweiten Art von Block zu starten.

Für noch mehr Leistung, möchten Sie vielleicht verschiedene Blocktypen hinzufügen und eine größere Größe des Blocks. Wenn die Blöcke immer die gleiche Größe und groß genug sind, können Sie sogar weniger Bits für die Zeiger verwenden. Mit einer Blockgröße von 16 Zeigern, haben Sie bereits einen Byte frei in einem 32-Bit-Adressraum. Werfen Sie einen Blick auf die Judy Baum Dokumentation für interessante Blocktypen. Grundsätzlich Sie Code und Engineering-Zeit für einen Raum (und Laufzeit) trade-off

hinzufügen

Sie wollen wahrscheinlich mit einem 256 breiten direkten Radix für die ersten vier Zeichen starten. Das bietet einen angemessenen Raum / Zeit-Kompromiss. In dieser Implementierung bekommen Sie sehr viel weniger Speicher-Overhead als mit einem einfachen Trie; es ist etwa dreimal kleiner (ich habe nicht gemessen). O (n) ist kein Problem, wenn die Konstante niedrig genug ist, wie Sie bemerkt, wenn sie mit dem O (n log n) quicksort verglichen wird.

Haben Sie Interesse verdoppelt im Umgang? Mit kurzen Sequenzen, es gehen zu sein. die Blöcke Anpassung zählt zu handhaben ist schwierig, aber es kann sehr platzeffizient sein.

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top