Je n'ai jamais vu une sorte de radix en place, et de la nature de la sorte de radix je doute qu'il soit beaucoup plus rapide qu'une sorte de out of place aussi longtemps que le tableau temporaire s'inscrit dans la mémoire.

raison:

le tri fait une lecture linéaire sur le tableau d'entrée, mais toutes les Écritures seront presque aléatoires. À partir d'un certain N vers le haut, cela se résume à une mise en cache de miss per write. Cette erreur de cache est ce qui ralentit votre algorithme. Si il est en place ou pas ne changera pas cet effet.

je sais que cela ne répondra pas directement à votre question, mais si le tri est un goulot d'étranglement, vous pouvez vouloir jeter un oeil à près du tri algorithmes comme un étape de prétraitement (le wiki-page sur le soft-heap peut vous aider à commencer).

qui pourrait donner un très bon coup de pouce à la localisation de cache. Un tri radix de Livre de texte hors de la place sera alors plus performante. Le les Écritures seront encore presque aléatoires mais au moins elles se regrouperont autour des mêmes morceaux de mémoire et en tant que tel augmentent le taux de succès de cache.

je n'ai aucune idée si cela fonctionne dans la pratique.

Btw: si vous avez affaire à des chaînes D'ADN seulement: vous pouvez compresser un char en deux bits et emballer vos données beaucoup. Cela réduira l'exigence de mémoire du facteur quatre par rapport à une représentation naïve. L'adressage devient plus complexe, mais L'ALU de de toute façon, votre CPU a beaucoup de temps à passer pendant toutes les pannes de cache.

basé sur le code de dsimcha, j'ai implémenté une version plus générique qui s'intègre bien dans le cadre que nous utilisons (SeqAn). En fait, transférer le code était très simple. Ce n'est qu'après que j'ai trouvé qu'il y avait qui sont en fait des publications concernant ce même sujet. Ce qui est génial, c'est qu'ils disent la même chose que vous. Un article d'Andersson et Nilsson sur Implementing Radixsort vaut certainement la peine d'être lu. S'il vous arrive de connaissez l'Allemand, assurez-vous également de lire la thèse de diplôme de David Weese où il met en œuvre un indice générique de substrat. La plupart de la thèse est consacrée à une analyse détaillée du coût de construction de l'indice, compte tenu secondaire de la mémoire et des fichiers extrêmement volumineux. Les résultats de son travail ont effectivement été mis en œuvre dans SeqAn, mais pas dans les parties où j'en avais besoin.

juste pour le plaisir, voici le code que j'ai écrit (Je ne pense pas que quiconque n'utilise pas SeqAn n'aura aucune utilité). Notez qu'il ne considère toujours pas radixes plus grand 4. Je m'attends à ce que cela ait un impact énorme sur la performance, mais malheureusement, je n'ai tout simplement pas le temps de le mettre en œuvre maintenant.

le code exécute plus de deux fois plus vite que Introsort pour les chaînes courtes. Même le point de rupture est à une longueur d'environ 12-13. Le type de chaîne (par exemple si elle a 4, 5, ou 16 valeurs différentes) est comparativement sans importance. Tri > 6 000 000 ADN les lectures du chromosome 2 du génome humain prennent un peu plus de 2 secondes sur mon PC. Juste pour info, c'est fast ! Surtout si l'on considère que je n'utilise pas SIMD ou toute autre accélération matérielle. De plus, valgrind me montre que le goulot d'étranglement principal est operator new dans les assignations de chaîne. Il est appelé environ 65 millions de fois-dix fois pour chaque corde! C'est un indice que swap pourrait être optimisé pour ces chaînes: au lieu de faire des copies, il on pourrait échanger tous les personnages. Je n'ai pas essayé, mais je suis convaincu que ça ferait une sacrée différence. Et, juste pour le dire à nouveau, au cas où quelqu'un n'écoutait pas: la taille radix n'a presque aucune influence sur l'exécution – ce qui signifie que je devrais certainement essayer de mettre en œuvre la suggestion faite par FryGuy, Stephan et EvilTeach.

Ah oui, au fait: la localité de cache est un facteur notable : à partir de chaînes de 1m, la durée d'exécution n'augmente plus de façon linéaire. Cependant, cela pourrait être corrigé assez facilement: j'utilise le tri d'insertion pour les petits sous – ensembles (<= 20 chaînes) - au lieu de mergesort comme suggéré par le hacker aléatoire. Apparemment, cela fonctionne encore mieux que mergesort pour de telles petites listes (voir le premier papier i lié).

namespace seqan {

template <typename It, typename F, typename T>
inline void prescan(It front, It back, F op, T const& id) {
    using namespace std;
    if (front == back) return;
    typename iterator_traits<It>::value_type accu = *front;
    *front++ = id;
    for (; front != back; ++front) {
        swap(*front, accu);
        accu = op(accu, *front);
    }
}

template <typename TIter, typename TSize, unsigned int RADIX>
inline void radix_permute(TIter front, TIter back, TSize (& bounds)[RADIX], TSize base) {
    for (TIter i = front; i != back; ++i)
        ++bounds[static_cast<unsigned int>((*i)[base])];

    TSize fronts[RADIX];

    std::copy(bounds, bounds + RADIX, fronts);
    prescan(fronts, fronts + RADIX, std::plus<TSize>(), 0);
    std::transform(bounds, bounds + RADIX, fronts, bounds, plus<TSize>());

    TSize active_base = 0;

    for (TIter i = front; i != back; ) {
        if (active_base == RADIX - 1)
            return;
        while (fronts[active_base] >= bounds[active_base])
            if (++active_base == RADIX - 1)
                return;
        TSize current_base = static_cast<unsigned int>((*i)[base]);
        if (current_base <= active_base)
            ++i;
        else
            std::iter_swap(i, front + fronts[current_base]);
        ++fronts[current_base];
    }
}

template <typename TIter, typename TSize>
inline void insertion_sort(TIter front, TIter back, TSize base) {
    typedef typename Value<TIter>::Type T;
    struct {
        TSize base, len;
        bool operator ()(T const& a, T const& b) {
            for (TSize i = base; i < len; ++i)
                if (a[i] < b[i]) return true;
                else if (a[i] > b[i]) return false;
            return false;
        }
    } cmp = { base, length(*front) }; // No closures yet. :-(

    for (TIter i = front + 1; i != back; ++i) {
        T value = *i;
        TIter j = i;
        for ( ; j != front && cmp(value, *(j - 1)); --j)
            *j = *(j - 1);
        if (j != i)
            *j = value;
    }
}

template <typename TIter, typename TSize, unsigned int RADIX>
inline void radix(TIter top, TIter front, TIter back, TSize base, TSize (& parent_bounds)[RADIX], TSize next) {
    if (back - front > 20) {
        TSize bounds[RADIX] = { 0 };
        radix_permute(front, back, bounds, base);

        // Sort current bucket recursively by suffix.
        if (base < length(*front) - 1)
            radix(front, front, front + bounds[0], base + 1, bounds, static_cast<TSize>(0));
    }
    else if (back - front > 1)
        insertion_sort(front, back, base);

    // Sort next buckets on same level recursively.
    if (next == RADIX - 1) return;
    radix(top, top + parent_bounds[next], top + parent_bounds[next + 1], base, parent_bounds, next + 1);
}

template <typename TIter>
inline void radix_sort(TIter front, TIter back) {
    typedef typename Container<TIter>::Type TStringSet;
    typedef typename Value<TStringSet>::Type TString;
    typedef typename Value<TString>::Type TChar;
    typedef typename Size<TStringSet>::Type TSize;

    TSize const RADIX = ValueSize<TChar>::VALUE;
    TSize bounds[RADIX];

    radix(front, front, back, static_cast<TSize>(0), bounds, RADIX - 1);
}

} // namespace seqan

vous pouvez certainement supprimer les exigences de mémoire en encodant la séquence en bits. Vous regardez des permutations donc, pour la longueur 2, avec "ACGT" c'est 16 états, ou 4 bits. Pour la longueur 3, c'est 64 États, qui peuvent être encodés en 6 bits. Donc ça ressemble à 2 bits pour chaque lettre dans la séquence, ou environ 32 bits pour 16 caractères comme vous l'avez dit.

S'il existe un moyen de réduire le nombre de "mots" valides, une nouvelle compression peut être possible.

donc pour les séquences de longueur 3, on pourrait créer 64 seaux, peut-être de taille uint32, ou uint64. Les initialiser à zéro. Iterate à travers votre très grande liste de 3 séquences de char, et de les encoder comme ci-dessus. Utilisez ceci comme un indice, et incrémentez ce seau.

Répétez ceci jusqu'à ce que toutes vos séquences aient été traitées.

ensuite, régénérez votre liste.

itérer à travers les 64 seaux dans l'ordre, pour le nombre trouvé dans ce seau, de générer de nombreux cas de la séquence représentée par ce seau.

lorsque toutes les seaux ont été réitéré, vous avez votre tableau trié.

une séquence de 4, ajoute 2 bits, donc il y aurait 256 seaux. Une séquence de 5, ajoute 2 bits, donc il y aurait 1024 seaux.

à un moment donné, le nombre de seaux approchera vos limites. Si vous lisez les séquences d'un fichier, au lieu de les garder dans de mémoire, plus de mémoire disponible pour les seaux.

je pense que ce serait plus rapide que de faire le tri in situ car les seaux sont susceptibles de s'adapter dans votre ensemble de travail.

voici un hack qui montre la technique

#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, '"151900920"', 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;
}

si votre ensemble de données est si grand, alors je pense qu'une approche de tampon basée sur le disque serait la meilleure:

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

j'expérimenterais aussi le groupement en un plus grand nombre de seaux, par exemple, si votre chaîne était:

GATTACA

le premier appel MSB retournerait le seau pour GATT (256 seaux au total), de cette façon vous feriez moins de branches du tampon basé sur le disque. Cela peut ou ne peut pas améliorer la performance, donc expérimenter avec elle.

je vais prendre un risque et vous suggérer de passer à un tas/ tas-port implémentation. Cette suggestion s'accompagne de quelques hypothèses:

1. Vous contrôlez la lecture des données
2. Vous pouvez faire quelque chose de significatif avec les données triées dès que vous 'démarrer' obtenir un tri.

la beauté du tri tas/tas est que vous pouvez construire le tas pendant que vous lisez les données, et vous pouvez commencer à obtenir des résultats dès que vous avez construit le tas.

reculons. Si vous êtes si chanceux que vous pouvez lire les données de manière asynchrone (c'est - à-dire que vous pouvez poster une sorte de demande de lecture et être notifié quand certaines données sont prêtes), et puis vous pouvez construire une partie du tas en attendant que le prochain morceau de données à venir-même à partir du disque. Souvent, cette approche peut enfouir la plupart du coût de la moitié de votre Tri derrière le temps passé à obtenir le données.

une fois les données lues, le premier élément est déjà disponible. Selon l'endroit où vous envoyez les données, cela peut être génial. Si vous l'envoyez à un autre lecteur asynchrone, ou à un autre modèle d'événement parallèle, ou à L'interface utilisateur, vous pouvez envoyer des morceaux et des morceaux au fur et à mesure.

cela dit-si vous n'avez aucun contrôle sur la façon dont les données sont lues, et qu'elles sont lues synchrones, et que vous n'avez aucune utilité pour les données triées jusqu'à ce qu'elles soient entièrement écrites - ignorer tout cela. : (

voir les articles de Wikipedia:

• aéroport
• tas Binaire

du point de vue de la Performance, vous pourriez vouloir regarder un plus général des algorithmes de tri de comparaison de chaîne.

Actuellement vous finissez par toucher chaque élément de chaque corde, mais vous pouvez faire mieux!

en particulier, un burst sort est un très bon ajustement pour cette affaire. En prime, puisque burstsort est basé sur des essais, il fonctionne ridiculement bien pour les petites tailles de l'alphabet utilisées dans L'ADN / RNA, puisque vous n'avez pas besoin de construire une sorte d'un noeud de recherche ternaire, d'un hachage ou d'un autre schéma de compression de noeud trie dans l'implémentation trie. Les essais peuvent être utiles pour votre but final de type suffix-array.

une implémentation correcte de burstsort est disponible sur Source forge à http://sourceforge.net/projects/burstsort / -mais il n'est pas en place.

à des fins De comparaison, Le C-burstsort de mise en œuvre de couverts à http://www.cs.mu.oz.au/~rsinha / papers / SinhaRingZobel-2006.pdf benchmarks 4-5x plus rapide que quicksort et Radix tries pour certaines charges de travail typiques.

vous pourriez essayer d'utiliser un trie . Le tri des données est simplement itérer à travers l'ensemble de données et l'insérer; la structure est naturellement triée, et vous pouvez penser qu'elle est similaire à un arbre B (sauf qu'au lieu de faire des comparaisons, vous toujours utilisez des pointeurs indirects).

le comportement de mise en cache favorisera tous les noeuds internes, donc vous ne pourrez probablement pas améliorer cela; mais vous pouvez jouer avec le facteur de ramification de votre trie aussi bien (assurez-vous que chaque noeud s'insère dans une ligne de cache unique, attribuez des noeuds de trie similaires à un tas, comme un tableau contigu qui représente un transversal d'ordre de niveau). Puisque les essais sont aussi des structures numériques (O (k) insert/find/delete pour les éléments de longueur k), vous devriez avoir des performances compétitives à un tri radix.

je burstsort un pique-bits de la représentation des chaînes de caractères. Burstsort est censé avoir une meilleure localisation que radix tries, ce qui réduit l'utilisation de l'espace supplémentaire par burst tries à la place des essais classiques. Le papier original a des mesures.

vous voulez jeter un oeil à Large-scale Genome Sequence Processing par les Drs Kasahara et Morishita.

composées des quatre lettres nucléotidiques A, C, G et T peuvent être spécialement encodées en nombres entiers pour un traitement plus rapide de much . Radix sort est parmi les nombreux algorithmes discutés dans le livre; vous devriez être en mesure d'adapter la réponse acceptée à cette question et de voir une grande amélioration de la performance.

le tri Radix sans espace supplémentaire " est un papier qui traite de votre problème.

Radix-Sort n'est pas conscient du cache et n'est pas l'algorithme de tri le plus rapide pour les grands ensembles. Vous pouvez regarder:

• ti7qsort . ti7qsort est le tri le plus rapide pour les entiers (peut être utilisé pour les chaînes de petite taille).
• Inline QSORT
• Chaîne de tri

vous pouvez également utiliser la compression et encoder chaque lettre de votre ADN en 2 bits avant de stocker dans le tableau de tri.

dsimcha's MSB radix sort semble agréable, mais Nils se rapproche au cœur du problème avec l'observation que la localisation de cache est ce qui vous tue à de grandes tailles de problèmes.

je suggère une approche très simple:

1. estimation Empirique la plus grande taille m , pour lequel une base de tri est efficace.
2. lire des blocs de m éléments à la fois, les Trier par radix, et les écrire (à un mémoire tampon si vous avez assez de mémoire, mais sinon de fichier), jusqu'à épuisement de votre entrée.
3. Mergesort l'résultant triés blocs.

Mergesort est l'algorithme de tri le plus convivial pour le cache que je connaisse:" Lire l'élément suivant du tableau A ou B, puis écrire un élément dans le tampon de sortie."Il fonctionne efficacement sur les lecteurs de bande . Il faut de l'espace 2n pour trier n items, mais mon pari est que la localisation de cache grandement améliorée que vous verrez rendra cela sans importance -- et si vous utilisiez une sorte de radix non-en-place, vous aviez besoin de cet espace supplémentaire de toute façon.

veuillez noter enfin que mergesort peut être implémenté sans récursion, et en fait le faire de cette façon permet de clarifier le vrai motif linéaire d'accès à la mémoire.

on dirait que vous avez résolu le problème, mais pour info, il semble qu'une version d'un tri radix en place est le"tri Drapeau Américain". C'est décrit ici: Engineering Radix Sort . L'idée générale est de faire 2 passes sur chaque caractère - d'abord compter combien de chaque Vous avez, de sorte que vous pouvez subdiviser le tableau d'entrée en bins. Puis passez à nouveau, en échangeant chaque élément dans la corbeille appropriée. Maintenant triez récursivement chaque corbeille sur la suivante position du personnage.

pensez D'abord au codage de votre problème. Débarrassez-vous des chaînes, remplacez-les par une représentation binaire. Utilisez le premier octet pour indiquer length+encoding. Alternativement, utilisez une représentation de longueur fixe à une limite de quatre octets. Alors le tri radix devient beaucoup plus facile. Pour un tri radix, la chose la plus importante est de ne pas avoir de manipulation d'exception au point chaud de la boucle intérieure.

OK, j'ai réfléchi un peu plus au problème des 4-naires. Vous voulez une solution comme un Judy arbre pour cela. La solution suivante peut gérer des chaînes de longueur variable; pour une longueur fixe, il suffit de supprimer les bits de longueur, ce qui rend en fait plus facile.

attribuent des blocs de 16 pointeurs. La partie la moins significative des pointeurs peut être réutilisée, car vos blocs seront toujours alignés. Vous pourriez vouloir un allocator de stockage spécial pour cela (fractionnant le grand stockage en blocs plus petits). Il existe différents types de blocs:

• Encodage avec 7 longueur de bits de chaînes de longueur variable. Au fur et à mesure qu'ils se remplissent, vous les remplacez par:
• position code les deux caractères suivants, vous avez 16 pointeurs vers les blocs suivants, se terminant par:
• Bitmap encodage des trois derniers caractères d'une chaîne.

pour chaque type de bloc, vous devez stocker des informations différentes dans les LSBs. Comme vous l'avez chaînes de longueur variable, vous devez stockez la fin de chaîne aussi, et le dernier type de bloc ne peut être utilisé que pour les cordes les plus longues. Les 7 bits de longueur doivent être remplacés par moins que vous obtenez plus profondément dans la structure.

cela vous fournit un stockage raisonnablement rapide et très efficace en mémoire des chaînes triées. Il se comportera un peu comme un trie . Pour que cela fonctionne, assurez-vous de construire suffisamment de tests unitaires. Vous voulez couvrir toutes les transitions de blocs. Vous voulez commencer avec seulement le deuxième type de bloc.

Pour encore plus de performances, vous pouvez ajouter différents types de blocs et une plus grande taille de bloc. Si les blocs sont toujours de la même taille et assez grand, vous pouvez utiliser encore moins de bits pour les pointeurs. Avec une taille de bloc de 16 pointeurs, vous avez déjà un octet de libre dans un espace d'adressage 32 bits. Jetez un oeil à la documentation de Judy tree pour les types de blocs intéressants. En gros, vous ajoutez du code et du temps d'ingénierie pour un compromis espace (et temps d'exécution)

vous voulez probablement commencer avec un radix direct de 256 wide pour les quatre premiers caractères. Cela fournit un compromis espace/temps décent. Dans cette implémentation, vous obtenez beaucoup moins de mémoire aérienne qu'avec un simple test; il est environ trois fois plus petit (je n'ai pas mesuré). O(n) n'est pas un problème si la constante est assez basse, comme vous l'avez remarqué en comparant avec O (N log n) quicksort.

êtes-vous intéressé par les doubles? Avec manches les séquences, il va y avoir. Adapter les blocs à la poignée de compte est délicat, mais il peut être très efficace.