Ceci est un long texte. S'il vous plaît, supportez-moi. En résumé, la question est: existe-t-il un algorithme de tri Radix sur place ?

Préliminaire

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

Pour le moment, j'utilise std::sortqui utilise introsort dans toutes les implémentations courantes de la STL . Cela fonctionne assez bien. Cependant, je suis convaincu que le tri Radix correspond parfaitement à mon problème et devrait fonctionner beaucoup 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 rapide). Cependant, le temps d'exécution se dégrade de façon catastrophique lorsque la taille du problème augmente ( N > 5 000 000).

La raison est évidente: le tri radix nécessite de copier toutes les données (plus d'une fois dans mon implémentation naïve, en fait). Cela signifie que j'ai mis ~ 4 Gio dans ma mémoire principale, ce qui tue évidemment les performances. Même si ce n'était pas le cas, je ne peux pas me permettre d'utiliser autant de mémoire car les tailles de problème deviennent encore plus importantes.

Cas d'utilisation

Idéalement, cet algorithme devrait fonctionner avec n'importe quelle longueur de chaîne entre 2 et 100, pour l'ADN ainsi que pour l'ADN5 (qui permet un caractère générique supplémentaire "N"), ou même l'ADN avec des codes d'ambiguïté IUPAC (résultant en 16 valeurs distinctes). Cependant, je me rends compte que tous ces cas ne peuvent pas être couverts, donc je suis satisfait de toute amélioration de vitesse que j'obtiens. Le code peut décider dynamiquement vers quel algorithme envoyer.

Recherche

Malheureusement, l'article Wikipédia sur le tri radix est inutile. La section sur une variante sur place est une poubelle complète. La section NIST-DADS sur le tri radix est pratiquement inexistante. Il existe un article à consonance prometteuse appelé Efficient Adaptive In-Place Radix Sorting qui décrit l'algorithme «MSL». Malheureusement, ce document est également décevant.

En particulier, il y a les choses suivantes.

Premièrement, l'algorithme contient plusieurs erreurs et laisse beaucoup inexpliqué. En particulier, il ne détaille pas l'appel de récursivité (je suppose simplement qu'il incrémente ou réduit un certain pointeur pour calculer les valeurs de décalage et de masque actuelles). De plus, il utilise les fonctions dest_groupet dest_addresssans donner de définitions. Je ne vois pas comment les implémenter efficacement (c'est-à-dire dans O (1); au moins, ce dest_addressn'est pas trivial).

Enfin et surtout, l'algorithme atteint la place en échangeant des indices de tableau avec des éléments à l'intérieur du tableau d'entrée. Cela ne fonctionne évidemment que sur des tableaux numériques. Je dois l'utiliser sur des cordes. Bien sûr, je pourrais juste taper un typage fort et continuer en supposant que la mémoire tolérera que je stocke un index où 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 des entiers 32 bits). Cela ne fait que 16 caractères (ignorons pour l'instant que 16> log (5 000 000)).

Un autre article de l'un des auteurs ne donne aucune description précise, mais il donne l'exécution de MSL comme sous-linéaire, ce qui est complètement faux.

Pour récapituler : Y a-t-il un espoir de trouver une implémentation de référence de travail ou au moins un bon pseudocode / description d'un type de radix en place qui fonctionne sur les chaînes d'ADN?

Eh bien, voici une implémentation simple d'un tri radix MSD pour l'ADN. Il est écrit en D parce que c'est la langue que j'utilise le plus et est donc moins susceptible de faire des erreurs idiotes, mais il pourrait facilement être traduit dans une autre langue. Il est en place mais nécessite des 2 * seq.lengthpassages dans 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, c'est un peu spécifique à l'ADN, plutôt que d'être général, mais cela devrait être rapide.

Éditer:

Je suis curieux de savoir si ce code fonctionne réellement, alors je l'ai testé / débogué en attendant que mon propre code bioinformatique s'exécute. La version ci-dessus est actuellement testée et fonctionne. Pour 10 millions de séquences de 5 bases chacune, c'est environ 3 fois plus rapide qu'un introsort optimisé.

Je n'ai jamais vu de tri radix sur place, et de par la nature du tri radix, je doute qu'il soit beaucoup plus rapide qu'un tri hors place tant que le tableau temporaire tient en 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 à un échec de cache par écriture. Cette erreur de cache est ce qui ralentit votre algorithme. S'il est en place ou non, cela 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 examiner les algorithmes de tri proches comme une étape de prétraitement (la page wiki sur le tas logiciel peut vous aider à démarrer).

Cela pourrait donner un très bon coup de pouce à la localisation du cache. Un tri radix hors-texte des manuels sera alors plus performant. Les écritures seront toujours presque aléatoires mais au moins, elles se regrouperont autour des mêmes morceaux de mémoire et augmenteront ainsi le taux d'accès au cache.

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

Btw: Si vous traitez uniquement avec des chaînes d'ADN: vous pouvez compresser un caractère en deux bits et emballer vos données beaucoup. Cela réduira les besoins en mémoire du facteur quatre sur 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 ratés de cache.

Vous pouvez certainement supprimer les besoins en mémoire en encodant la séquence en bits. Vous regardez les 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. 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 des séquences de longueur 3, on pourrait créer 64 compartiments, peut-être de taille uint32 ou uint64. Initialisez-les à zéro. Parcourez votre très très grande 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 afin, pour le nombre trouvé dans ce compartiment, de générer autant d'instances de la séquence représentée par ce compartiment.
lorsque tous les compartiments ont été itérés, vous disposez de 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 compartiments.

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

Je pense que ce serait plus rapide que de faire le tri in situ car les godets sont susceptibles de s'adapter à 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, '\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;
}

Si votre ensemble de données est si volumineux, je pense qu'une approche de tampon basée 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

J'expérimenterais également le regroupement en un plus grand nombre de compartiments, par exemple, si votre chaîne était:

GATTACA

le premier appel MSB renverrait le compartiment pour GATT (256 compartiments au total), de cette façon vous faites moins de branches du tampon basé sur le disque. Cela peut ou non améliorer les performances, alors essayez-les.

Je vais sortir sur un membre et vous suggère de passer à un tas / 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 trier.

La beauté du tas / tri en tas est que vous pouvez créer le tas pendant que vous lisez les données, et vous pouvez commencer à obtenir des résultats au moment où vous avez construit le tas.

Revenons en arrière. Si vous êtes si chanceux que vous pouvez lire les données de manière asynchrone (c'est-à-dire, vous pouvez publier une sorte de demande de lecture et être averti lorsque certaines données sont prêtes), puis vous pouvez créer un morceau du tas pendant que vous attendez le prochain bloc de données à venir - même à partir du disque. Souvent, cette approche peut enterrer la majeure partie du coût de la moitié de votre tri derrière le temps passé à obtenir 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 parfait. Si vous l'envoyez à un autre lecteur asynchrone, ou à un modèle d'événement ou d'interface utilisateur parallèle, 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 de manière synchrone, et que vous n'avez aucune utilité pour les données triées jusqu'à ce qu'elles soient entièrement écrites - ignorez tout cela. :(

Voir les articles Wikipedia:

• Heapsort
• Tas binaire

" Tri Radix sans espace supplémentaire " est un document qui résout votre problème.

En termes de performances, vous souhaiterez peut-être examiner des algorithmes de tri de comparaison de chaînes plus généraux.

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

En particulier, un tri en rafale convient très bien à ce cas. En prime, puisque burstsort est basé sur des essais, cela fonctionne ridiculement bien pour les petites tailles d'alphabet utilisées dans l'ADN / ARN, car vous n'avez pas besoin de construire une sorte de nœud de recherche ternaire, de hachage ou autre schéma de compression de nœud de tri dans le mise en œuvre de trois. Les essais peuvent également être utiles pour votre objectif final de type tableau de suffixes.

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

À des fins de comparaison, l'implémentation de C-burstsort a couvert à http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf les tests de référence 4 à 5 fois plus rapidement que les types quicksort et radix pour certaines charges de travail typiques.

Vous aurez envie de jeter un œil au traitement des séquences de génomes à grande échelle par les Drs. Kasahara et Morishita.

Les chaînes composées des quatre lettres nucléotidiques A, C, G et T peuvent être spécialement codées en nombres entiers pour un traitement beaucoup plus rapide. Le tri Radix fait partie des 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 des performances.

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

Le comportement de mise en cache favorisera tous les nœuds internes, donc vous n'améliorerez probablement pas 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 de tri similaires à un tas, comme un tableau contigu qui représente une traversée d'ordre de niveau). Étant donné que les essais sont également des structures numériques (O (k) insert / find / delete pour les éléments de longueur k), vous devriez avoir des performances compétitives par rapport à un tri radix.

Je voudrais éclater une représentation compacte des chaînes. Burstsort est censé avoir une bien meilleure localité que les sortes de radix, ce qui réduit l'utilisation d'espace supplémentaire avec des essais en rafale à la place des essais classiques. Le papier d'origine a des dimensions.

Radix-Sort n'est pas sensible au 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 fixe).
• Inline QSORT
• Tri des chaînes

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

Le tri radix MSB de dsimcha a l'air bien, mais Nils se rapproche du cœur du problème avec l'observation que la localité du cache est ce qui vous tue à des problèmes de grande taille.

Je propose une approche très simple:

1. Estimez empiriquement la plus grande taille mpour laquelle un tri radix est efficace.
2. Lisez des blocs d' méléments à la fois, triez-les par radix et écrivez-les (dans un tampon de mémoire si vous avez suffisamment de mémoire, mais sinon à classer), jusqu'à épuisement de votre entrée.
3. Fusionne les blocs triés résultants.

Mergesort est l'algorithme de tri le plus convivial pour le cache que je connaisse: "Lisez l'élément suivant dans le tableau A ou B, puis écrivez un élément dans le tampon de sortie." Il fonctionne efficacement sur les lecteurs de bande . Cela nécessite de l' 2nespace pour trier les néléments, mais je parie que la localité de cache beaucoup améliorée que vous verrez rendra cela sans importance - et si vous utilisiez un tri radix non en place, vous aviez besoin de cet espace supplémentaire de toute façon.

Veuillez noter enfin que le mergesort peut être implémenté sans récursivité, et en fait le faire de cette façon rend clair le véritable modèle d'accès à la mémoire linéaire.

Il semble que vous ayez résolu le problème, mais pour mémoire, il semble qu'une version d'un tri Radix en place réalisable soit le "American Flag Sort". Il est décrit ici: Engineering Radix Sort . L'idée générale est de faire 2 passes sur chaque caractère - comptez d'abord combien vous en avez, afin de pouvoir subdiviser le tableau d'entrée en cases. Puis recommencez, en échangeant chaque élément dans le bon bac. Triez maintenant récursivement chaque casier sur la position de caractère suivante.

Tout d'abord, pensez 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 + le codage. Vous pouvez également utiliser une représentation de longueur fixe à une limite de quatre octets. Ensuite, le tri radix devient beaucoup plus facile. Pour un tri radix, le plus important est de ne pas avoir de gestion d'exception au point chaud de la boucle interne.

OK, j'ai réfléchi un peu plus au problème des 4 naires. Vous voulez une solution comme un arbre Judy 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 facilite la tâche.

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 cela (diviser le grand stockage 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 encode les deux caractères suivants, vous avez 16 pointeurs vers 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 différentes informations 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 à mesure que vous approfondissez la structure.

Cela vous offre 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 générer suffisamment de tests unitaires. Vous voulez une couverture de toutes les transitions de bloc. Vous souhaitez commencer avec uniquement le deuxième type de bloc.

Pour encore plus de performances, vous souhaiterez peut-être ajouter différents types de blocs et une plus grande taille de bloc. 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 disposez déjà d'un octet libre dans un espace d'adressage 32 bits. Jetez un œil à la documentation de l'arborescence Judy pour les types de blocs intéressants. Fondamentalement, vous ajoutez du code et du temps d'ingénierie pour un compromis d'espace (et d'exécution)

Vous voudrez probablement commencer avec un radix direct de 256 larges pour les quatre premiers caractères. Cela fournit un compromis espace / temps décent. Dans cette implémentation, vous obtenez beaucoup moins de surcharge de mémoire qu'avec un simple trie; 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 O (n log n).

Êtes-vous intéressé à gérer des doubles? Avec de courtes séquences, il va y en avoir. L'adaptation des blocs pour gérer les nombres est délicate, mais elle peut être très économe en espace.