Ordinamento Radix sul posto

StackOverflow https://stackoverflow.com/questions/463105

Domanda

Questo è un testo lungo. Per favore abbi pazienza. In breve, la domanda è: Esiste un algoritmo di ordinamento radix sul posto praticabile ?

preliminare

Ho un numero enorme di stringhe di piccole dimensioni a lunghezza fissa che utilizzano solo le lettere & # 8220; A & # 8221 ;, & # 8220; C & # 8221 ;, & # 8220; G & # 8221; e & # 8220; T & # 8221; (sì, l'hai indovinato: DNA ) che desidero ordinare.

Al momento, utilizzo std :: sort che utilizza introsort in tutte le implementazioni comuni di STL . Funziona abbastanza bene. Tuttavia, sono convinto che ordinamento radix si adatta perfettamente al mio problema e dovrebbe funzionare molto meglio in pratica.

Dettagli

Ho testato questo presupposto con un'implementazione molto ingenua e per input relativamente piccoli (dell'ordine di 10.000) questo era vero (beh, almeno più del doppio della velocità). Tuttavia, il tempo di esecuzione si riduce in modo abissale quando la dimensione del problema aumenta ( N > 5.000.000).

Il motivo è ovvio: l'ordinamento radix richiede la copia di tutti i dati (più di una volta nella mia ingenua implementazione, in realtà). Ciò significa che ho inserito ~ 4 GiB nella mia memoria principale, il che ovviamente uccide le prestazioni. Anche se non fosse così, non posso permettermi di usare così tanta memoria poiché le dimensioni del problema in realtà diventano ancora più grandi.

Usa casi

Idealmente, questo algoritmo dovrebbe funzionare con qualsiasi lunghezza di stringa compresa tra 2 e 100, sia per il DNA che per il DNA5 (che consente un carattere jolly aggiuntivo & # 8220; N & # 8221;) o persino il DNA con IUPAC codici di ambiguità (risultanti in 16 valori distinti). Tuttavia, mi rendo conto che tutti questi casi non possono essere coperti, quindi sono contento di qualsiasi miglioramento della velocità che ottengo. Il codice può decidere dinamicamente a quale algoritmo inviare.

Ricerca

Sfortunatamente, l'articolo di Wikipedia sull'ordinamento radix è inutile. La sezione relativa a una variante sul posto è spazzatura completa. La sezione NIST-DADS sull'ordinamento radix è quasi inesistente. C'è un documento dal suono promettente chiamato Efficient Adaptive In-Place Radix Sorting che descrive l'algoritmo & # 8220; MSL & # 8221 ;. Sfortunatamente, anche questo documento è deludente.

In particolare, ci sono le seguenti cose.

Innanzitutto, l'algoritmo contiene diversi errori e lascia molto inspiegabile. In particolare, non dettaglia la chiamata di ricorsione (presumo semplicemente che aumenti o riduca alcuni puntatori per calcolare i valori di spostamento e maschera correnti). Inoltre, utilizza le funzioni dest_group e dest_address senza fornire definizioni. Non riesco a vedere come implementarli in modo efficiente (vale a dire in O (1); almeno dest_address non è banale).

Ultimo ma non meno importante, l'algoritmo raggiunge la posizione sul posto scambiando gli indici dell'array con elementi all'interno dell'array di input. Questo ovviamente funziona solo su array numerici. Ho bisogno di usarlo sulle stringhe. Certo, potrei semplicemente rovinare la digitazione forte e andare avanti supponendo che la memoria tollererà la mia memorizzazione di un indice in cui non appartiene. Ma questo funziona solo finché riesco a comprimere le mie stringhe in 32 bit di memoria (assumendo numeri interi a 32 bit). Sono solo 16 caratteri (ignoriamo per il momento che 1

È stato utile?

Soluzione

Bene, ecco una semplice implementazione di un ordinamento radix MSD per il DNA. È scritto in D perché è la lingua che uso di più e quindi è meno probabile che commetta errori sciocchi, ma potrebbe facilmente essere tradotto in un'altra lingua. È installato ma richiede 2 * seq.length attraverso l'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);
   }
}

Ovviamente, questo è un po 'specifico per il DNA, al contrario di essere generale, ma dovrebbe essere veloce.

Modifica:

Mi sono incuriosito se questo codice funziona davvero, quindi l'ho testato / debug mentre aspettavo che il mio codice bioinformatico fosse eseguito. La versione sopra ora è attualmente testata e funziona. Per 10 milioni di sequenze di 5 basi ciascuna, è circa 3 volte più veloce di un introsort ottimizzato.

Altri suggerimenti

Non ho mai visto un ordinamento radix sul posto, e dalla natura dell'ordinamento radix dubito che sia molto più veloce di un ordinamento fuori posto fintanto che l'array temporaneo si adatta alla memoria.

Motivo:

L'ordinamento esegue una lettura lineare sull'array di input, ma tutte le scritture saranno quasi casuali. Da un certo N in poi questo si riduce a una cache mancata per scrittura. Questa mancata cache è ciò che rallenta il tuo algoritmo. Se è a posto o no non cambierà questo effetto.

So che questo non risponderà direttamente alla tua domanda, ma se l'ordinamento è un collo di bottiglia potresti voler dare un'occhiata agli algoritmi near sorting come fase di preelaborazione ( la pagina wiki sul soft-heap potrebbe farti iniziare).

Ciò potrebbe dare un ottimo impulso alla localizzazione della cache. Un ordinamento radix fuori posto da manuale funzionerà quindi meglio. Le scritture saranno ancora quasi casuali ma almeno si raggrupperanno attorno agli stessi blocchi di memoria e di conseguenza aumenteranno il rapporto di hit della cache.

Non ho idea se funzionerà in pratica però.

Btw: se hai a che fare solo con stringhe di DNA: puoi comprimere un carattere in due bit e comprimere abbastanza i tuoi dati. Ciò ridurrà il fabbisogno di memoria del fattore quattro rispetto a una rappresentazione naiiva. L'indirizzamento diventa più complesso, ma l'ALU della tua CPU ha comunque molto tempo da spendere durante tutti i mancati cache.

Puoi certamente eliminare i requisiti di memoria codificando la sequenza in bit. Stai osservando le permutazioni quindi, per la lunghezza 2, con "ACGT" sono 16 stati o 4 bit. Per la lunghezza 3, sono 64 stati, che possono essere codificati in 6 bit. Quindi sembrano 2 bit per ogni lettera nella sequenza, o circa 32 bit per 16 caratteri come hai detto.

Se esiste un modo per ridurre il numero di "parole" valide, potrebbe essere possibile un'ulteriore compressione.

Quindi, per sequenze di lunghezza 3, si potrebbero creare 64 bucket, magari di dimensione uint32 o uint64. Inizializzali a zero. Scorri il tuo elenco molto ampio di 3 sequenze di caratteri e codificali come sopra. Usa questo come un pedice e incrementa quel bucket.
Ripeti fino a quando tutte le tue sequenze non sono state elaborate.

Successivamente, rigenera la tua lista.

Scorrere i 64 bucket in ordine, per il conteggio trovato in quel bucket, genera così tante istanze della sequenza rappresentata da quel bucket.
quando tutti i bucket sono stati ripetuti, hai il tuo array ordinato.

Una sequenza di 4, aggiunge 2 bit, quindi ci sarebbero 256 bucket. Una sequenza di 5, aggiunge 2 bit, quindi ci sarebbero 1024 bucket.

Ad un certo punto il numero di bucket si avvicinerà ai tuoi limiti. Se leggi le sequenze da un file, invece di tenerle in memoria, sarebbe disponibile più memoria per i bucket.

Penso che questo sarebbe più veloce che fare l'ordinamento in situ poiché è probabile che i secchi si adattino al tuo set di lavoro.

Ecco un hack che mostra la tecnica 

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

Se il tuo set di dati è così grande, penso che un approccio buffer basato su disco sarebbe il migliore: 

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

Sperimenterei anche il raggruppamento in un numero maggiore di bucket, ad esempio, se la stringa fosse: 

GATTACA

la prima chiamata MSB restituisce il bucket per GATT (256 bucket totali), in questo modo si creano meno rami del buffer basato su disco. Questo può o meno migliorare le prestazioni, quindi sperimentale.

Per quanto riguarda le prestazioni, potresti voler esaminare algoritmi di ordinamento di confronto delle stringhe più generali.

Attualmente finisci per toccare ogni elemento di ogni stringa, ma puoi fare di meglio!

In particolare, un burst sort è molto adatto per questo caso. Come bonus, dal momento che burstsort si basa sui tentativi, funziona ridicolmente bene per le piccole dimensioni dell'alfabeto utilizzate nel DNA / RNA, poiché non è necessario creare alcun tipo di nodo di ricerca ternario, hash o altro schema di compressione del nodo trie nel trie implementazione. I tentativi possono essere utili anche per il tuo obiettivo finale simile a un array di suffissi.

Un'implementazione decente di burstsort per scopi generici è disponibile su forge di origine all'indirizzo http://sourceforge.net/projects / burstsort / - ma non è sul posto.

A fini di confronto, l'implementazione di C-burstsort è stata trattata in http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf benchmark 4-5 volte più veloce di quicksort e ordinamenti radix per alcuni carichi di lavoro tipici.

Ti consigliamo di dare un'occhiata a Elaborazione sequenziale del genoma su larga scala di Drs. Kasahara e Morishita.

Le stringhe composte dalle quattro lettere nucleotidiche A, C, G e T possono essere appositamente codificate in numeri interi per un'elaborazione molto più rapida. L'ordinamento Radix è tra i molti algoritmi discussi nel libro; dovresti essere in grado di adattare la risposta accettata a questa domanda e vedere un notevole miglioramento delle prestazioni.

" Ordinamento Radix senza spazio extra " è un documento che affronta il tuo problema.

Potresti provare a utilizzare un trie . L'ordinamento dei dati sta semplicemente ripetendo il set di dati e inserendolo; la struttura è naturalmente ordinata e puoi considerarla simile a un B-Tree (tranne che invece di fare confronti, sempre usi le indicazioni indirette del puntatore).

Il comportamento della cache favorirà tutti i nodi interni, quindi probabilmente non migliorerai su quello; ma puoi anche giocherellare con il fattore di ramificazione del tuo trie (assicurati che ogni nodo si adatti a una singola riga della cache, alloca nodi trie simili a un heap, come un array contiguo che rappresenta un attraversamento di ordine di livello). Poiché i tentativi sono anche strutture digitali (O (k) inserisci / trova / elimina per elementi di lunghezza k), dovresti avere prestazioni competitive per un ordinamento radix.

Vorrei burstsort una rappresentazione a bit delle stringhe. Si dice che Burstsort abbia una localizzazione molto migliore rispetto a quella dei radix, mantenendo basso lo spazio aggiuntivo con i tentativi di scoppio al posto dei tentativi classici. Il documento originale ha misure.

Radix-Sort non è consapevole della cache e non è l'algoritmo di ordinamento più veloce per set di grandi dimensioni. Puoi guardare:

Puoi anche usare la compressione e codificare ogni lettera del tuo DNA in 2 bit prima di archiviarli nella matrice di ordinamento.

L'ordinamento radix MSB di dsimcha sembra carino, ma Nils si avvicina al nocciolo del problema con l'osservazione che la localizzazione della cache è ciò che ti sta uccidendo con problemi di grandi dimensioni.

Suggerisco un approccio molto semplice:

  1. Stima empiricamente la dimensione m di dimensioni maggiori per cui un ordinamento radix è efficiente.
  2. Leggi i blocchi di m alla volta, radixali e scrivili (su un buffer di memoria se hai memoria sufficiente, ma altrimenti su file), fino a esaurire l'input.
  3. Fusione i blocchi ordinati risultanti.

Mergesort è l'algoritmo di ordinamento più adatto alla cache di cui sono a conoscenza: " Leggi l'elemento successivo dall'array A o B, quindi scrivi un elemento nel buffer di output. " Funziona in modo efficiente su unità nastro . Richiede lo spazio 2n per ordinare gli elementi n , ma la mia scommessa è che la localizzazione della cache molto migliorata che vedrai renderà poco importante - e se stessi usando un ordinamento radix non sul posto, avevi comunque bisogno di quello spazio extra.

Si noti infine che il mergesort può essere implementato senza ricorsione, e infatti farlo in questo modo chiarisce il vero modello di accesso alla memoria lineare.

Sembra che tu abbia risolto il problema, ma per la cronaca, sembra che una versione di un ordinamento radix sul posto funzionante sia il "quotazione della bandiera americana". È descritto qui: Engineering Radix Sort . L'idea generale è di fare 2 passaggi su ogni personaggio - prima conta quanti di ciascuno di essi hai, quindi puoi suddividere l'array di input in bin. Quindi ripassare, scambiando ogni elemento nel cestino corretto. Ora ordina in modo ricorsivo ogni cestino nella posizione successiva del personaggio.

Innanzitutto, pensa alla codifica del tuo problema. Sbarazzarsi delle stringhe, sostituirle con una rappresentazione binaria. Utilizzare il primo byte per indicare lunghezza + codifica. In alternativa, utilizzare una rappresentazione a lunghezza fissa con un limite di quattro byte. Quindi l'ordinamento radix diventa molto più semplice. Per un ordinamento radix, la cosa più importante è non avere la gestione delle eccezioni nel punto caldo del ciclo interno.

OK, ho pensato un po 'di più al problema 4-nary. Per questo, vuoi una soluzione come un Judy tree . La soluzione successiva può gestire stringhe di lunghezza variabile; per lunghezza fissa basta rimuovere i bit di lunghezza, che in realtà lo rendono più facile.

Alloca blocchi di 16 puntatori. Il bit meno significativo dei puntatori può essere riutilizzato, poiché i blocchi saranno sempre allineati. Potrebbe essere necessario un allocatore di memoria speciale per esso (suddividere la memoria di grandi dimensioni in blocchi più piccoli). Esistono diversi tipi di blocchi:

  • Codifica con 7 bit di lunghezza di stringhe di lunghezza variabile. Mentre si riempiono, li sostituisci con:
  • Posizione codifica i due caratteri successivi, hai 16 puntatori ai blocchi successivi, che terminano con:
  • Codifica bitmap degli ultimi tre caratteri di una stringa.

Per ogni tipo di blocco, è necessario memorizzare informazioni diverse negli LSB. Dato che hai stringhe di lunghezza variabile, devi archiviare anche end-of-string e l'ultimo tipo di blocco può essere utilizzato solo per le stringhe più lunghe. I 7 bit di lunghezza dovrebbero essere sostituiti da meno man mano che approfondisci la struttura.

Ciò fornisce una memorizzazione ragionevolmente veloce e molto efficiente della memoria delle stringhe ordinate. Si comporterà in qualche modo come un trie . Per farlo funzionare, assicurati di costruire abbastanza unit test. Vuoi copertura di tutte le transizioni di blocco. Vuoi iniziare solo con il secondo tipo di blocco.

Per prestazioni ancora maggiori, potresti voler aggiungere diversi tipi di blocco e una dimensione del blocco più grande. Se i blocchi hanno sempre le stesse dimensioni e sono abbastanza grandi, puoi usare ancora meno bit per i puntatori. Con una dimensione di blocco di 16 puntatori, hai già un byte libero in uno spazio di indirizzi a 32 bit. Dai un'occhiata alla documentazione dell'albero di Judy per i tipi di blocchi interessanti. Fondamentalmente, aggiungi il codice e il tempo di progettazione per un compromesso spaziale (e di runtime)

Probabilmente vuoi iniziare con una radix diretta larga 256 per i primi quattro caratteri. Ciò fornisce un discreto compromesso spazio / tempo. In questa implementazione, si ottiene molto meno sovraccarico di memoria rispetto a un semplice trie; è circa tre volte più piccolo (non ho misurato). O (n) non è un problema se la costante è abbastanza bassa, come hai notato durante il confronto con il quicksort O (n log n).

Sei interessato a gestire i doppi? Con brevi sequenze, ci saranno. Adattare i blocchi per gestire i conteggi è complicato, ma può essere molto efficiente in termini di spazio.

Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top