To jest długi tekst. Proszę o wyrozumiałość. Sprowadzone pytanie brzmi: czy istnieje praktyczny algorytm sortowania radix w miejscu ?

Wstępny

Mam ogromną liczbę małych ciągów o stałej długości, które używają tylko liter „A”, „C”, „G” i „T” (tak, zgadłeś: DNA ), które chcę posortować.

W tej chwili używam, std::sortktóry wykorzystuje introsort we wszystkich popularnych implementacjach STL . To działa całkiem dobrze. Jestem jednak przekonany, że sortowanie radix idealnie pasuje do mojego zestawu problemów i powinno działać znacznie lepiej w praktyce.

Detale

Przetestowałem to założenie z bardzo naiwną implementacją i przy stosunkowo niewielkich nakładach (rzędu 10 000) było to prawdą (cóż, przynajmniej dwa razy szybciej). Jednak środowisko wykonawcze obniża się gwałtownie, gdy rozmiar problemu staje się większy ( N > 5 000 000).

Powód jest oczywisty: sortowanie radix wymaga skopiowania całych danych (tak naprawdę więcej niż raz w mojej naiwnej implementacji). Oznacza to, że umieściłem ~ 4 GiB w mojej głównej pamięci, co oczywiście zabija wydajność. Nawet jeśli nie, nie mogę sobie pozwolić na użycie tak dużej ilości pamięci, ponieważ rozmiary problemów stają się jeszcze większe.

Przypadków użycia

Idealnie, ten algorytm powinien działać z dowolną długością łańcucha od 2 do 100, zarówno dla DNA, jak i DNA5 (co pozwala na dodatkowy znak wieloznaczny „N”), a nawet DNA z kodami niejednoznaczności IUPAC (co daje 16 różnych wartości). Zdaję sobie jednak sprawę, że nie można uwzględnić wszystkich tych przypadków, więc cieszę się z każdej poprawy prędkości, jaką otrzymuję. Kod może dynamicznie decydować, do którego algorytmu wysłać.

Badania

Niestety artykuł Wikipedii na temat sortowania radix jest bezużyteczny. Część dotycząca wariantu na miejscu to kompletne śmieci. Sekcja NIST-DADS na temat sortowania radix jest prawie nieistniejąca. Istnieje obiecująco brzmiący artykuł o nazwie Efficient Adaptive In-Place Radix Sorting, który opisuje algorytm „MSL”. Niestety, ten artykuł również rozczarowuje.

W szczególności są następujące rzeczy.

Po pierwsze, algorytm zawiera kilka błędów i pozostawia wiele niewyjaśnionych. W szczególności nie wyszczególnia wywołania rekurencyjnego (po prostu zakładam, że zwiększa lub zmniejsza wskaźnik, aby obliczyć bieżące wartości przesunięcia i maski). Korzysta także z funkcji dest_groupi dest_addressnie podaje definicji. Nie widzę, jak efektywnie je wdrożyć (to znaczy w O (1); przynajmniej dest_addressnie jest to trywialne).

Na koniec algorytm osiąga miejsce w miejscu, zamieniając indeksy tablic na elementy wewnątrz tablicy wejściowej. To oczywiście działa tylko na tablice numeryczne. Muszę go używać na ciągach. Oczywiście mógłbym po prostu mocno wkręcić i pisać dalej, zakładając, że pamięć będzie tolerować przechowywanie indeksu, do którego on nie należy. Ale to działa tylko tak długo, jak długo mogę wycisnąć moje ciągi do 32 bitów pamięci (zakładając 32-bitowe liczby całkowite). To tylko 16 znaków (zignorujmy na razie, że 16> log (5 000 000)).

Kolejny artykuł jednego z autorów nie zawiera żadnego dokładnego opisu, ale podaje środowisko wykonawcze MSL jako sublinearne, co jest całkowicie błędne.

Podsumowując : Czy jest jakaś nadzieja na znalezienie działającej implementacji referencyjnej lub przynajmniej dobrego pseudokodu / opisu działającego na miejscu sortowania radix, który działa na łańcuchach DNA?

Oto prosta implementacja sortowania MSD dla DNA. Jest napisany w D, ponieważ jest to język, którego używam najczęściej i dlatego najmniej prawdopodobne jest, że popełnię głupie błędy, ale można go łatwo przetłumaczyć na inny język. Jest na miejscu, ale wymaga 2 * seq.lengthprzejścia przez tablicę.

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

Oczywiście jest to swoiste dla DNA, a nie ogólne, ale powinno być szybkie.

Edytować:

Ciekawe, czy ten kod rzeczywiście działa, więc przetestowałem go / debugowałem, czekając na uruchomienie mojego kodu bioinformatycznego. Powyższa wersja jest obecnie testowana i działa. Dla 10 milionów sekwencji po 5 zasad każda jest około 3 razy szybsza niż zoptymalizowana introsort.

Nigdy nie widziałem sortowania radix na miejscu, a z natury sortowania radix wątpię, aby było ono znacznie szybsze niż sortowanie nie na miejscu, o ile tablica tymczasowa mieści się w pamięci.

Powód:

Sortowanie dokonuje liniowego odczytu tablicy wejściowej, ale wszystkie zapisy będą prawie losowe. Od pewnego N w górę sprowadza się to do braku pamięci podręcznej na zapis. Ten brak pamięci podręcznej spowalnia Twój algorytm. Jeśli jest na miejscu, czy nie, nie zmieni tego efektu.

Wiem, że to nie odpowie bezpośrednio na twoje pytanie, ale jeśli sortowanie jest wąskim gardłem, możesz przyjrzeć się algorytmom blisko sortowania jako krokowi wstępnego przetwarzania (strona wiki na miękkim stosie może zacząć).

To może dać bardzo ładny wzrost lokalizacji pamięci podręcznej. Sortowanie według podręcznika w miejscu poza miejscem będzie wtedy działać lepiej. Zapisy nadal będą prawie losowe, ale przynajmniej skupią się wokół tych samych fragmentów pamięci i jako takie zwiększą współczynnik trafień w pamięci podręcznej.

Nie mam jednak pojęcia, czy to zadziała w praktyce.

Btw: Jeśli masz do czynienia tylko z ciągami DNA: możesz skompresować znak do dwóch bitów i spakować swoje dane całkiem sporo. To zmniejszy zapotrzebowanie na pamięć czterokrotnie w stosunku do naiwnej reprezentacji. Adresowanie staje się bardziej złożone, ale ALU twojego procesora i tak ma dużo czasu do spędzenia podczas wszystkich braków pamięci podręcznej.

Z pewnością możesz zmniejszyć wymagania dotyczące pamięci, kodując sekwencję w bitach. Patrzysz na permutacje, więc dla długości 2 z „ACGT”, który ma 16 stanów lub 4 bity. Dla długości 3 jest to 64 stany, które można zakodować w 6 bitach. Wygląda więc na 2 bity na każdą literę w sekwencji lub około 32 bity na 16 znaków, jak powiedziałeś.

Jeśli istnieje sposób na zmniejszenie liczby prawidłowych „słów”, możliwa jest dalsza kompresja.

Tak więc dla sekwencji o długości 3 można utworzyć 64 wiadra, może mieć rozmiar uint32 lub uint64. Zainicjuj je do zera. Iteruj po swojej bardzo dużej liście 3 sekwencji znaków i koduj je jak wyżej. Użyj tego jako indeksu dolnego i zwiększaj ten segment.
Powtarzaj to do momentu przetworzenia wszystkich sekwencji.

Następnie ponownie wygeneruj listę.

Iteruj po 64 segmentach w kolejności, aby uzyskać liczbę znalezioną w tym segmencie, wygeneruj tyle wystąpień sekwencji reprezentowanych przez to segment.
gdy wszystkie segmenty zostały iterowane, masz posortowaną tablicę.

Sekwencja 4, dodaje 2 bity, więc będzie 256 wiader. Sekwencja 5, dodaje 2 bity, więc będzie 1024 wiader.

W pewnym momencie liczba wiader zbliży się do twoich limitów. Jeśli odczytasz sekwencje z pliku, zamiast przechowywać je w pamięci, dostępna będzie większa pamięć dla segmentów.

Myślę, że byłoby to szybsze niż robienie tego na miejscu, ponieważ wiadra prawdopodobnie mieszczą się w twoim zestawie roboczym.

Oto hack, który pokazuje technikę

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

Jeśli twój zestaw danych jest tak duży, pomyślałbym, że najlepszym rozwiązaniem byłoby zastosowanie bufora dyskowego:

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

Eksperymentowałbym również grupowanie w większą liczbę segmentów, na przykład, jeśli Twój ciąg był:

GATTACA

pierwsze wywołanie MSB zwróci segment dla GATT (256 całkowitych segmentów), w ten sposób utworzysz mniej gałęzi bufora opartego na dysku. To może, ale nie musi, poprawić wydajność, więc eksperymentuj z tym.

Mam zamiar wyjść na kończynę i zasugerować przejście na implementację heap / heapsort . Ta sugestia zawiera pewne założenia:

1. Kontrolujesz odczyt danych
2. Możesz zrobić coś znaczącego z posortowanymi danymi, gdy tylko zaczniesz je sortować.

Piękno sterty / sortowania sterty polega na tym, że można zbudować stertę podczas odczytywania danych, a wyniki można zacząć od momentu zbudowania sterty.

Cofnijmy się. Jeśli masz tyle szczęścia, że ​​możesz odczytać dane asynchronicznie (tzn. Możesz opublikować jakieś żądanie odczytu i otrzymać powiadomienie, gdy niektóre dane będą gotowe), a następnie możesz zbudować część sterty, czekając na następna porcja danych, która ma wejść - nawet z dysku. Często takie podejście może pogrzebać większość kosztów połowy sortowania w stosunku do czasu poświęconego na uzyskanie danych.

Po odczytaniu danych pierwszy element jest już dostępny. W zależności od miejsca przesyłania danych może to być świetne. Jeśli wysyłasz go do innego asynchronicznego czytnika lub innego równoległego modelu „zdarzenia” lub interfejsu użytkownika, możesz wysyłać porcje i porcje w trakcie pracy.

To powiedziawszy - jeśli nie masz kontroli nad tym, jak dane są odczytywane, a dane są odczytywane synchronicznie, a posortowane dane nie są używane, dopóki nie zostaną całkowicie zapisane - zignoruj ​​to wszystko. :(

Zobacz artykuły w Wikipedii:

• Heapsort
• Kupa binarna

„ Sortowanie Radix bez dodatkowej przestrzeni ” to artykuł rozwiązujący Twój problem.

Pod względem wydajności warto przyjrzeć się bardziej ogólnym algorytmom sortowania porównań ciągów.

Obecnie kończysz dotykając każdego elementu każdego sznurka, ale możesz zrobić to lepiej!

W szczególności rodzaj serii jest bardzo dobrze dopasowany do tego przypadku. Jako bonus, ponieważ burstsort opiera się na próbach, działa absurdalnie dobrze dla małych rozmiarów alfabetu używanych w DNA / RNA, ponieważ nie trzeba budować żadnego trójskładnikowego węzła wyszukiwania, skrótu lub innego schematu kompresji węzła trie w wdrożenie trie. Próby te mogą być również przydatne do ostatecznego celu podobnego do tablicy przyrostków.

Przyzwoite ogólne zastosowanie burstsort jest dostępne na stronie source forge pod adresem http://sourceforge.net/projects/burstsort/ - ale nie ma go na miejscu.

Dla celów porównawczych implementacja C-burstsort opisana na stronie http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf testy porównawcze 4-5 razy szybsze niż sortowanie Quicksort i Radix dla niektórych typowych obciążeń.

Będziesz chciał przyjrzeć się przetwarzaniu sekwencji genomu na dużą skalę przez Drs. Kasahara i Morishita.

Ciągi złożone z czterech liter nukleotydowych A, C, G i T mogą być specjalnie kodowane w liczbach całkowitych w celu znacznie szybszego przetwarzania. Sortowanie Radix jest jednym z wielu algorytmów omówionych w książce; powinieneś być w stanie dostosować przyjętą odpowiedź do tego pytania i zobaczyć znaczną poprawę wydajności.

Możesz spróbować użyć trie . Sortowanie danych polega na iteracji zestawu danych i wstawianiu go; struktura jest naturalnie posortowana i można ją traktować jako podobną do B-drzewa (z wyjątkiem tego, że zamiast dokonywania porównań, zawsze używa się pośrednich wskaźników).

Zachowanie w pamięci podręcznej sprzyja wszystkim wewnętrznym węzłom, więc prawdopodobnie nie poprawisz tego; ale możesz również manipulować współczynnikiem rozgałęzienia swojego trie (upewnij się, że każdy węzeł mieści się w jednej linii pamięci podręcznej, alokuj węzły trie podobne do sterty, jako ciągłą tablicę reprezentującą przechodzenie przez kolejność poziomów). Ponieważ próby są również strukturami cyfrowymi (O (k) wstaw / znajdź / usuń dla elementów o długości k), powinieneś mieć konkurencyjną wydajność do sortowania radix.

Chciałbym burstsort reprezentację pakowane-bitowej strun. Twierdzi się, że Burstsort ma znacznie lepszą lokalizację niż rodzaje radix, dzięki czemu dodatkowe użycie przestrzeni jest mniejsze dzięki próbom seryjnym zamiast próbom klasycznym. Oryginalny papier ma wymiary.

Sortowanie Radix nie obsługuje pamięci podręcznej i nie jest najszybszym algorytmem sortowania dla dużych zestawów. Możesz spojrzeć na:

• ti7qsort . ti7qsort to najszybszy sposób sortowania liczb całkowitych (może być używany do łańcuchów o małych, stałych rozmiarach).
• Inline QSORT
• Sortowanie ciągów

Możesz także użyć kompresji i zakodować każdą literę swojego DNA na 2 bity przed zapisaniem w tablicy sortowania.

Sortowanie MSB dsimcha wygląda ładnie, ale Nils zbliża się do sedna problemu, obserwując, że lokalizacja pamięci podręcznej zabija cię przy dużych rozmiarach problemu.

Proponuję bardzo proste podejście:

1. Oszacuj empirycznie największy rozmiar, mdla którego sortowanie radix jest wydajne.
2. Odczytuj bloki melementów na raz, sortuj je radix i zapisuj (do bufora pamięci, jeśli masz wystarczającą ilość pamięci, ale w innym przypadku do pliku), aż do wyczerpania danych wejściowych.
3. Połącz uzyskane posortowane bloki.

Mergesort jest najbardziej przyjaznym dla pamięci podręcznej algorytmem sortowania, jaki znam: „Odczytaj następny element z tablicy A lub B, a następnie zapisz element do bufora wyjściowego”. Działa wydajnie na napędach taśmowych . Nie wymaga 2nmiejsca do sortowania nprzedmiotów, ale założę się, że znacznie ulepszona lokalizacja pamięci podręcznej, którą zobaczysz, sprawi, że nie będzie to ważne - a jeśli używałeś sortowania radix w miejscu, i tak potrzebujesz dodatkowej przestrzeni.

Na koniec zauważ, że scalesort może zostać zaimplementowany bez rekurencji, aw rzeczywistości wykonanie tego w ten sposób wyjaśnia prawdziwy wzorzec dostępu do pamięci liniowej.

Wygląda na to, że rozwiązałeś problem, ale dla przypomnienia wydaje się, że jedną z wersji praktycznego sortowania radix w miejscu jest „American Flag Sort”. Jest to opisane tutaj: Inżynieria Sortowanie Radix . Ogólna idea polega na wykonaniu 2 przejść dla każdego znaku - najpierw policz, ile masz każdego z nich, abyś mógł podzielić tablicę wejściową na przedziały. Następnie przejdź ponownie, zamieniając każdy element do właściwego pojemnika. Teraz rekurencyjnie posortuj każdy pojemnik na następnej pozycji postaci.

Najpierw pomyśl o kodowaniu swojego problemu. Pozbądź się ciągów, zamień je na reprezentację binarną. Użyj pierwszego bajtu, aby wskazać długość + kodowanie. Alternatywnie użyj stałej reprezentacji długości na granicy czterech bajtów. Następnie sortowanie radix staje się znacznie łatwiejsze. W przypadku rodzaju radix najważniejszą rzeczą jest brak obsługi wyjątków w gorącym punkcie wewnętrznej pętli.

OK, myślałem trochę więcej o czwartym problemie. Potrzebujesz rozwiązania takiego jak drzewo Judy . Następne rozwiązanie może obsługiwać łańcuchy o zmiennej długości; dla stałej długości wystarczy usunąć bity długości, co faktycznie ułatwia.

Przydziel bloki po 16 wskaźników. Najmniej znaczącą część wskaźników można ponownie wykorzystać, ponieważ bloki zawsze będą wyrównane. Możesz potrzebować specjalnego alokatora pamięci (dzielenie dużej pamięci na mniejsze bloki). Istnieje wiele różnych rodzajów bloków:

• Kodowanie za pomocą 7 bitów ciągów o zmiennej długości. Gdy się zapełniają, zastępujesz je:
• Pozycja koduje kolejne dwa znaki, masz 16 wskaźników do następnych bloków, kończących się na:
• Kodowanie bitmapowe ostatnich trzech znaków ciągu.

Dla każdego rodzaju bloku musisz przechowywać różne informacje w LSB. Ponieważ masz ciągi o zmiennej długości, musisz również przechowywać koniec łańcucha, a ostatniego rodzaju bloku można używać tylko dla najdłuższych ciągów. 7 bitów długości powinno zostać zastąpionych przez mniej, gdy wejdziesz głębiej w strukturę.

Zapewnia to dość szybkie i bardzo wydajne pamięciowe sortowanie ciągów znaków. Zachowuje się trochę jak trie . Aby to zadziałało, należy zbudować wystarczającą liczbę testów jednostkowych. Chcesz objąć wszystkie przejścia blokowe. Chcesz zacząć od drugiego rodzaju bloku.

Aby uzyskać jeszcze większą wydajność, możesz chcieć dodać różne typy bloków i większy rozmiar bloku. Jeśli bloki są zawsze tego samego rozmiaru i wystarczająco duże, możesz użyć jeszcze mniej bitów dla wskaźników. Przy rozmiarze bloku 16 wskaźników masz już bajt wolny w 32-bitowej przestrzeni adresowej. Przejrzyj dokumentację drzewa Judy, aby znaleźć interesujące typy bloków. Zasadniczo dodajesz kod i czas inżynierii dla kompromisu przestrzeni (i środowiska wykonawczego)

Prawdopodobnie chcesz zacząć od bezpośredniej podstawy o szerokości 256 dla pierwszych czterech znaków. To zapewnia przyzwoity kompromis czas / przestrzeń. W tej implementacji uzyskujesz znacznie mniej pamięci niż w przypadku zwykłej wersji próbnej; jest około trzy razy mniejszy (nie mierzyłem). O (n) nie stanowi problemu, jeśli stała jest wystarczająco niska, jak zauważyłeś podczas porównywania z szybkim sortowaniem O (n log n).

Czy jesteś zainteresowany obsługą podwójnych? Będą to krótkie sekwencje. Dostosowanie bloków do obsługi liczników jest trudne, ale może być bardzo wydajne pod względem miejsca.