Tri sur place sur place

Question

Ceci est un texte long. S'il vous plaît supporter avec moi. En résumé, la question qui se pose est la suivante: Existe-t-il un algorithme de tri de radix sur place utilisable ?

Préliminaire

J'ai un grand nombre de petites chaînes de longueur fixe qui utilisent uniquement les lettres & # 8220; A & # 8221 ;, C & # 8221 ;, & # 8220; G & # 8221; et & # 8220; T & # 8221; (oui, vous l'avez deviné: ADN ) que je veux trier.

Pour le moment, j'utilise std :: sort , qui utilise introsort dans toutes les implémentations courantes de la STL . Cela fonctionne assez bien. Cependant, je suis convaincu que le sorte de radix convient parfaitement à mon problème et qu'il devrait beaucoup mieux fonctionner mieux dans la pratique.

Détails

J'ai testé cette hypothèse avec une implémentation très naïve et pour des entrées relativement petites (de l'ordre de 10 000), cela était vrai (enfin, au moins plus de deux fois plus vite). Cependant, l'exécution se dégrade de manière catastrophique lorsque la taille du problème augmente ( N > 5 000 000).

La raison en est évidente: le tri de base nécessite la copie de l’ensemble des données (plus d’une fois dans mon implémentation naïve, en fait). Cela signifie que j'ai mis environ 4 Gio dans ma mémoire principale, ce qui tue évidemment les performances. Même si ce n’était pas le cas, je ne pourrais pas me permettre d’utiliser autant de mémoire, car la taille du problème devenait encore plus grande.

Cas d'utilisation

Idéalement, cet algorithme devrait fonctionner avec toute longueur de chaîne comprise entre 2 et 100, pour l'ADN aussi bien que pour l'ADN5 (ce qui permet d'ajouter un caractère générique supplémentaire), ou même un ADN avec IUPAC codes d'ambiguïté (donnant 16 valeurs distinctes). Cependant, je me rends compte que tous ces cas ne peuvent pas être couverts, je suis donc satisfait de toute amélioration de vitesse que je reçois Le code peut décider dynamiquement de l’algorithme à envoyer.

Recherche

Malheureusement, le article de Wikipedia sur le tri en radix est inutile. La section sur une variante en place est un déchet complet. La section NIST-DADS sur le tri de base est quasi inexistante. Il existe un article prometteur appelé Tri efficace sur radix sur place adaptable , qui décrit l'algorithme & # 8220; MSL & # 8221 ;. Malheureusement, cet article aussi est décevant.

En particulier, il y a les choses suivantes.

Premièrement, l’algorithme contient plusieurs erreurs et laisse beaucoup d’explications inexpliquées. En particulier, il ne décrit pas en détail l’appel récursif (je suppose simplement qu’il incrémente ou réduit un pointeur pour calculer les valeurs de décalage et de masque en cours). De plus, il utilise les fonctions groupe_dest et adresse_dest sans donner de définitions. Je ne vois pas comment les mettre en œuvre efficacement (c'est-à-dire dans O (1); au moins dest_address n'est pas trivial).

Dernier point, mais non le moindre, l’algorithme permet d’atteindre son emplacement en échangeant des index de tableau avec des éléments à l’intérieur du tableau en entrée. Cela ne fonctionne évidemment que sur les tableaux numériques. J'ai besoin de l'utiliser sur des chaînes. Bien sûr, je pourrais simplement taper du texte fort et continuer en supposant que la mémoire me permettra de stocker un index auquel il n’appartient pas. Mais cela ne fonctionne que tant que je peux compresser mes chaînes dans 32 bits de mémoire (en supposant que les entiers 32 bits). C’est seulement 16 caractères (ignorons pour le moment que 1

Answer 1

Eh bien, voici une implémentation simple d’un type de base MSD pour l’ADN. Il est écrit en D parce que c’est la langue que j’utilise le plus et que j’ai donc le moins de chances de faire des erreurs stupides, mais il pourrait facilement être traduit dans une autre langue. Il est en place mais nécessite 2 * seq.length passe par le tableau.

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

Évidemment, cela est plutôt spécifique à l'ADN, par opposition à la généralisation, mais cela devrait être rapide.

Modifier:

Je me demandais si ce code fonctionnait réellement. Je l’ai donc testé / débogué en attendant l’exécution de mon propre code de bioinformatique. La version ci-dessus est actuellement testée et fonctionne. Pour 10 millions de séquences de 5 bases chacune, il est environ 3 fois plus rapide qu'un introsort optimisé.

Answer 2

Je n'ai jamais vu une sorte de radix en place, et de par la nature de la sorte-radix, je doute qu'elle soit beaucoup plus rapide qu'un tri hors de propos tant que le tableau temporaire tient dans la mémoire.

Raison:

Le tri effectue une lecture linéaire sur le tableau d'entrée, mais toutes les écritures seront presque aléatoires. À partir d'un certain N, cela se résume à une cache manquante par écriture. Ce manque de cache est ce qui ralentit votre algorithme. Si c'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 voudrez peut-être jeter un coup d'œil aux algorithmes de tri proche en tant que étape de prétraitement ( la page wiki sur le tas souple peut vous aider à démarrer).

Cela pourrait donner un très bon coup de pouce à la localisation du cache. Un tri de base de texte non déplacé dans les manuels fonctionnera alors mieux. Les écritures seront toujours presque aléatoires, mais au moins elles se regrouperont autour des mêmes blocs de mémoire et augmenteront ainsi le taux d'accès au cache.

Je ne sais pas si cela fonctionne dans la pratique.

Btw: Si vous utilisez uniquement des chaînes d'ADN: vous pouvez compresser un caractère en deux bits et compresser vos données à volonté. Cela réduira les besoins en mémoire par un facteur quatre par rapport à une représentation naïve. L’adressage devient plus complexe, mais l’ALU de votre CPU a quand même beaucoup de temps à consacrer à tous les échecs en mémoire cache.

Answer 3

Vous pouvez certainement supprimer les besoins en mémoire en codant la séquence en bits. Vous recherchez des permutations donc, pour la longueur 2, avec " ACGT " c'est 16 états, ou 4 bits. Pour la longueur 3, cela représente 64 états, qui peuvent être codés en 6 bits. Cela ressemble donc à 2 bits pour chaque lettre de 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 compression supplémentaire peut être possible.

Ainsi, pour les séquences de longueur 3, vous pouvez créer 64 compartiments, de taille uint32 ou uint64. Initialisez les à zéro. Parcourez votre très très longue liste de 3 séquences de caractères et encodez-les comme ci-dessus. Utilisez-le comme indice et incrémentez ce compartiment.
Répétez cette opération jusqu'à ce que toutes vos séquences aient été traitées.

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

Parcourez les 64 compartiments dans l'ordre, pour le nombre trouvé dans ce compartiment, générez autant d'occurrences de la séquence représentée par ce compartiment.
lorsque tous les compartiments ont été itérés, vous avez votre tableau trié.

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

À un moment donné, le nombre de compartiments approchera de vos limites. Si vous lisez les séquences à partir d'un fichier, au lieu de les conserver en mémoire, davantage de mémoire sera disponible pour les compartiments.

Je pense que cela serait plus rapide que de faire le tri in situ car les seaux vont probablement s'intégrer dans votre environnement 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, '\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;
}

Answer 4

Si votre ensemble de données est si volumineux, alors je penserais qu’une approche basée sur un tampon basé sur 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

Je voudrais aussi expérimenter le regroupement dans un plus grand nombre de compartiments, par exemple, si votre chaîne était:

GATTACA

le premier appel MSB renverrait le compartiment pour le GATT (256 compartiments au total), ce qui réduirait le nombre de branches du tampon basé sur disque. Cela peut ou non améliorer les performances, alors testez-le.

Answer 5

Je vais vous dire quelque chose et vous suggérer de passer à un heap / heapsort mise en œuvre. 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 commencez à les faire trier.

La beauté du tas / tas-tri 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.

Faisons un pas en arrière. Si vous avez la chance de pouvoir lire les données de manière asynchrone (c’est-à-dire que vous pouvez envoyer une sorte de demande de lecture et être averti lorsque certaines données sont prêtes), vous pouvez ensuite créer une partie du tas pendant que vous attendez le message. prochain bloc de données à venir - même à partir du disque. Souvent, cette approche peut couvrir la majeure partie du coût de la moitié de votre tri après le temps passé à récupérer les 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 formidable. Si vous l'envoyez à un autre lecteur asynchrone, ou à un modèle "d'événement" parallèle, ou à une interface utilisateur, vous pouvez envoyer des fragments à la volée.

Cela dit - si vous ne contrôlez pas la façon dont les données sont lues, si elles sont lues de manière synchrone, et si vous n’utilisez plus les données triées tant qu’elles ne sont pas entièrement écrites - ignorez tout cela. : (

Voir les articles Wikipedia:

• Heapsort
• segment de mémoire binaire

Answer 6

En ce qui concerne les performances, vous pouvez envisager des algorithmes plus généraux de tri par comparaison de chaînes.

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

En particulier, un tri par rafales convient très bien dans ce cas. De plus, comme burstsort est basé sur des tentatives, il fonctionne parfaitement bien pour les alphabets de petite taille utilisés dans l’ADN / ARN, car il n’est pas nécessaire de créer un quelconque schéma de compression de nœud de recherche ternaire, de hachage ou un autre nœud. la mise en œuvre. Les essais peuvent également être utiles pour votre objectif final, semblable à un tableau de suffixes.

Une implémentation décente de burstsort à usage général est disponible sur source forge à l'adresse http://sourceforge.net/projects. / burstsort / - mais ce n'est pas sur place.

À des fins de comparaison, la mise en œuvre de C-burstsort est décrite à http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf benchmarks 4-5x plus rapidement que les tris rapides et tri-radix pour certaines charges de travail typiques.

Answer 7

Vous voudrez jeter un coup d’œil sur le Traitement de séquence de génome à grande échelle par les Drs. Kasahara et Morishita.

Les chaînes composées des quatre lettres de nucléotide A, C, G et T peuvent être spécialement codées en Integers pour un traitement beaucoup plus rapide. Le type Radix fait partie des nombreux algorithmes présentés dans le livre. vous devriez pouvoir adapter la réponse acceptée à cette question et constater une nette amélioration des performances.

Answer 8

" Tri de radix sans espace supplémentaire " est un document traitant de votre problème.

Answer 9

Vous pouvez essayer d'utiliser un trie . Le tri des données consiste simplement à parcourir l'ensemble de données et à l'insérer. la structure est naturellement triée et vous pouvez la considérer comme similaire à un arbre B (sauf qu'au lieu de faire des comparaisons, vous utilisez toujours des indirections de pointeur).

Le comportement de la mise en cache favorisera tous les nœuds internes, vous ne pourrez donc probablement pas améliorer cela; mais vous pouvez également jouer avec le facteur de branchement de votre trie (assurez-vous que chaque nœud tient dans une seule ligne de cache, allouez des nœuds semblables à un segment de mémoire, sous la forme d'un tableau contigu représentant une traversée ordre de niveau). Puisque les essais sont aussi des structures numériques (O (k) insert / find / delete pour des éléments de longueur k), vous devez avoir des performances compétitives comparées à un tri de base.

Answer 10

Je voudrais envoyer en rafale une représentation compacte des chaînes. Burstsort aurait une bien meilleure localisation que les tris radix, limitant l'utilisation d'espace supplémentaire par des essais en rafale au lieu d'essais classiques. Le papier d'origine contient des mesures.

Answer 11

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).
• QSORT en ligne
• Tri des chaînes

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

Answer 12

Le type de base MSB de dsimcha semble bien, mais Nils se rapproche du cœur du problème en observant que la localisation en cache est ce qui tue les problèmes de grande taille.

Je suggère une approche très simple:

1. Estimez empiriquement la plus grande taille m pour laquelle un tri de base est efficace.
2. Lisez des blocs d’éléments m à la fois, triez-les radicalement et écrivez-les (dans une mémoire tampon si vous avez assez de mémoire, mais sinon dans un fichier), jusqu’à épuisement de vos entrées.
3. Mergesort les blocs triés résultants.

Mergesort est l'algorithme de tri le plus convivial pour le cache que je connaisse: "Lisez le prochain élément du tableau A ou B, puis écrivez un élément dans le tampon de sortie." Il fonctionne efficacement sur les lecteurs de bande . Il faut un espace 2n pour trier les éléments n , mais je parie que la localité de cache bien améliorée que vous verrez rendra cela sans importance - et si vous utilisiez un tri de base 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écursivité. En procédant ainsi, vous définissez clairement le véritable modèle d'accès linéaire à la mémoire.

Answer 13

Il semble que vous ayez résolu le problème, mais aux fins du compte rendu, il apparaît qu'une version d'un type de tri de radix sur place utilisable est le "Tri par drapeau américain". Il est décrit ici: Tri par radier technique . L'idée générale est de faire 2 passages sur chaque caractère - commencez par compter le nombre de chacun, pour pouvoir subdiviser le tableau d'entrée en bacs. Recommencez ensuite, en échangeant chaque élément dans le bac approprié. Maintenant, triez chaque bac de manière récursive sur la position du caractère suivant.

Answer 14

Tout d’abord, réfléchissez 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 la longueur + codage. Vous pouvez également utiliser une représentation de longueur fixe à une limite de quatre octets. Ensuite, le tri de base devient beaucoup plus facile. Pour un tri à base, le plus important est de ne pas gérer les exceptions au point chaud de la boucle interne.

OK, j’ai réfléchi un peu plus au problème des 4-Nary. Vous souhaitez une solution comme un arbre Judy pour cela. La solution suivante peut gérer des chaînes de longueur variable; pour une longueur fixe, enlevez simplement les bits de longueur, ce qui facilite les choses.

Allouez des blocs de 16 pointeurs. Le bit le moins significatif des pointeurs peut être réutilisé, car vos blocs seront toujours alignés. Vous voudrez peut-être un allocateur de stockage spécial pour ce dernier (décomposer un stockage volumineux en blocs plus petits). Il existe différents types de blocs:

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

Pour chaque type de bloc, vous devez stocker des informations différentes dans les LSB. Comme vous avez des chaînes de longueur variable, vous devez également stocker la fin de chaîne, et le dernier type de bloc ne peut être utilisé que pour les chaînes les plus longues. Les 7 bits de longueur doivent être remplacés par moins au fur et à mesure que vous avancez dans la structure.

Vous disposez ainsi d’un stockage relativement 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 créer suffisamment de tests unitaires. Vous voulez une couverture de toutes les transitions de blocs. Vous voulez commencer avec le deuxième type de bloc uniquement.

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

Vous voudrez probablement commencer par une base directe large de 256 pour les quatre premiers caractères. Cela procure un compromis temps / espace décent. Dans cette implémentation, vous obtenez beaucoup moins de surcharge de mémoire 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 suffisamment basse, comme vous l’avez remarqué lors de la comparaison avec le tri rapide de O (n log n).

Êtes-vous intéressé par le traitement des doubles? Avec de courtes séquences, il y en aura. Il est difficile d’adapter les blocs pour gérer les comptes, mais cela peut être très économe en espace.