Tri par base sur place

Question

C'est un long texte. S'il vous plaît, soyez patient. En fin de compte, la question est: Existe-t-il un algorithme de tri de base en place fonctionnel ?

---

Préliminaire

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

En ce moment, j'utilise std::sort qui utilise l'algorithm introsort dans toutes les implémentations communes de la STL. Cela fonctionne plutôt bien. Cependant, je suis convaincu que le tri radix convient parfaitement à mon ensemble de problèmes et devrait fonctionner beaucoup mieux en 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), c'était vrai (environ deux fois plus rapide en réalité). Cependant, le temps d'exécution se dégrade de façon abyssale lorsque la taille du problème devient plus grande (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 réalité). Cela signifie que j'ai mis ~ 4 Gio dans ma mémoire principale ce qui tue évidemment les performances. Même si cela fonctionnait, je ne peux pas me permettre d'utiliser autant de mémoire car les tailles des problèmes deviennent en réalité encore plus grandes.

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 l'ADN5 (qui autorise un caractère de joker supplémentaire "N"), ou même l'ADN avec des codes d'ambiguïté IUPAC (ce qui donne 16 valeurs distinctes). Cependant, je comprends 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 de quel algorithme envoyer.

Recherche

Malheureusement, l'article Wikipedia sur le tri radix est inutile. La section sur une variante en place est complètement nulle. La section NIST-DADS sur le tri radix est quasi inexistante. Il y a un article prometteur appelé Tri radix en place efficace et adaptatif qui décrit l'algorithme "MSL". Malheureusement, cet article est également décevant.

En particulier, il y a les éléments suivants.

Tout d'abord, l'algorithme contient plusieurs erreurs et laisse beaucoup d'éléments non expliqués. En particulier, il ne détaille pas 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 actuelles). De plus, il utilise les fonctions dest_group et dest_address sans donner de définitions. Je ne vois pas comment implémenter efficacement ces fonctions (c'est-à-dire, en O(1) ; au moins dest_address n'est pas trivial).

Enfin, l'algorithme parvient à être en-place en échangeant les indices du tableau avec des éléments à l'intérieur du tableau d'entrée. Cela fonctionne évidemment uniquement sur les tableaux numériques. Je dois l'utiliser sur des chaînes. Bien sûr, je pourrais juste ignorer le typage fort et supposer que la mémoire tolérera que je stocke un indice là où il ne devrait pas être. Mais cela ne fonctionne que tant que je peux compresser mes chaînes en 32 bits de mémoire (en supposant des entiers sur 32 bits). Ce n'est que 16 caractères (ignorons pour l'instant que 16 > log(5 000 000)).

Un autre article par l'un des auteurs ne donne pas du tout de description précise, mais il indique un temps d'exécution de MSL inférieur à 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 fonctionnelle ou au moins un bon pseudo-code/description d'un tri radix in-place qui fonctionne sur des chaînes ADN ?

Answer 1

Eh bien, voici une implémentation simple d'un tri radix MSD pour l'ADN. Il est écrit en D car c'est le langage que j'utilise le plus et donc celui dans lequel je suis le moins susceptible de commettre des erreurs stupides, mais il pourrait facilement être traduit dans un autre langage. Il se fait sur place mais nécessite 2 * seq.length passages à travers 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);
   }
}

De toute évidence, cela est plutôt spécifique à l'ADN, par opposition à être général, mais cela devrait être rapide.

Éditer:

Je me suis demandé si ce code fonctionnait réellement, alors je l'ai testé/débogué tout en attendant que mon propre code de bioinformatique s'exécute. La version ci-dessus est en fait 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é.

Answer 2

Je n'ai jamais vu de tri radix in-place, et vu la nature du tri radix, je doute qu'il soit beaucoup plus rapide qu'un tri out-of-place tant que le tableau temporaire rentre en mémoire.

Raison :

Le tri effectue une lecture linéaire dans le tableau d'entrée, mais toutes les écritures seront presque aléatoires. À partir d'un certain N, cela revient à un manque de cache par écriture. Ce manque de cache est ce qui ralentit votre algorithme. S'il est in-place ou non 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 œil aux algorithmes de tri à proximité en tant qu'étape de prétraitement (la page wiki sur le soft-heap pourrait vous aider à démarrer).

Cela pourrait donner un très bel avantage en termes de localité du cache. Un tri radix out-of-place classique se comportera alors mieux. Les écritures seront toujours presque aléatoires mais au moins elles se regrouperont autour des mêmes parties de mémoire et augmenteront ainsi le taux de réussite du cache.

Je ne sais pas si cela fonctionne en pratique cependant.

Soit dit en passant : si vous ne travaillez qu'avec des chaînes d'ADN : vous pouvez compresser un caractère en deux bits et packer vos données considérablement. 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 beaucoup de temps à dépenser pendant tous les manques de cache de toute façon.

Answer 3

En s'appuyant sur le code de dsimcha, j'ai implémenté une version plus générique qui s'intègre bien dans le framework que nous utilisons (SeqAn). En réalité, porter le code était complètement simple. Ce n'est qu'ensuite que j'ai découvert qu'il existe effectivement des publications concernant ce même sujet. La bonne nouvelle est : elles disent essentiellement la même chose que vous. Un article d'Andersson et Nilsson sur La mise en œuvre de Radixsort vaut vraiment la peine d'être lu. Si vous connaissez l'allemand, assurez-vous également de lire la thèse de diplôme de David Weese où il met en œuvre un index de sous-chaîne générique. La plupart de la thèse est consacrée à une analyse détaillée du coût de la construction de l'index, en tenant compte de la mémoire secondaire et de fichiers extrêmement volumineux. Les résultats de son travail ont en fait été implémentés dans SeqAn, sauf dans les parties où j'en avais besoin.

Juste pour le plaisir, voici le code que j'ai écrit (je ne pense pas que quelqu'un qui n'utilise pas SeqAn en aura besoin). Remarquez qu'il ne prend toujours pas en compte les radices supérieurs à 4. Je m'attends à ce que cela ait un impact énorme sur les performances, mais malheureusement je n'ai tout simplement pas le temps en ce moment pour implémenter cela.

Le code s'exécute plus de deux fois plus vite que l'Introsort pour les chaînes courtes. Le point mort est à une longueur d'environ 12-13. Le type de chaîne (par exemple, qu'elle ait 4, 5 ou 16 valeurs différentes) est relativement peu important. Le tri de > 6 000 000 lectures d'ADN du chromosome 2 du génome humain prend un peu plus de 2 secondes sur mon PC. Juste pour préciser, c'est rapide ! Surtout en considérant que je n'utilise pas SIMD ou toute autre accélération matérielle. De plus, valgrind me montre que le principal goulot d'étranglement est operator new dans les affectations de chaînes. Il est appelé environ 65 000 000 fois - dix fois pour chaque chaîne ! C'est un signe révélateur que swap pourrait être optimisé pour ces chaînes : au lieu de faire des copies, il pourrait simplement échanger tous les caractères. Je n'ai pas essayé cela mais je suis convaincu que cela ferait une énorme différence. Et, juste pour le dire encore une fois, au cas où quelqu'un ne l'aurait pas entendu : la taille de la base presque aucun influence sur le temps d'exécution - ce qui signifie que je devrais définitivement essayer de mettre en œuvre la suggestion de FryGuy, Stephan et EvilTeach.

Ah oui, soit dit en passant : la localité du cache est un facteur notable : À partir de 1M de chaînes, le temps d'exécution n'augmente plus linéairement. Cependant, cela pourrait être corrigé assez facilement : j'utilise un tri par insertion pour de 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 listes de petite taille (voir le premier article que j'ai lié).

namespace seqan {

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

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

    TSize fronts[RADIX];

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

    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((*i)[base]);
        if (current_base <= active_base)
            ++i;
        else
            std::iter_swap(i, front + fronts[current_base]);
        ++fronts[current_base];
    }
}

template 
inline void insertion_sort(TIter front, TIter back, TSize base) {
    typedef typename Value::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 
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);

        // Trier de manière récursive le seau actuel par suffixe.
        if (base < length(*front) - 1)
            radix(front, front, front + bounds[0], base + 1, bounds, static_cast(0));
    }
    else if (back - front > 1)
        insertion_sort(front, back, base);

    // Trier de manière récursive les prochains seaux au même niveau.
    if (next == RADIX - 1) return;
    radix(top, top + parent_bounds[next], top + parent_bounds[next + 1], base, parent_bounds, next + 1);
}

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

    TSize const RADIX = ValueSize::VALUE;
    TSize bounds[RADIX];

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

} // namespace seqan

Answer 4

Vous pouvez certainement réduire les exigences de mémoire en encodant la séquence en bits. Vous regardez les permutations donc, pour une longueur de 2, avec "ACGT" c'est 16 états, soit 4 bits. Pour une longueur de 3, c'est 64 états, qui peuvent être encodés en 6 bits. Donc il semble qu'il y ait 2 bits pour chaque lettre dans la séquence, soit environ 32 bits pour 16 caractères comme vous l'avez dit.

S'il y a un moyen de réduire le nombre de 'mots' valides, une compression supplémentaire pourrait être possible.

Donc pour des séquences de longueur 3, on pourrait créer 64 seaux, peut-être de taille uint32, ou uint64. Initialisez-les à zéro. Parcourez votre très très grande liste de séquences de 3 caractères, et encodez-les comme indiqué ci-dessus. Utilisez ceci comme un indice, et incrémentez ce seau.
Répétez cela jusqu'à ce que toutes vos séquences aient été traitées.

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

Parcourez les 64 seaux dans l'ordre, pour le nombre trouvé dans ce seau, générez autant d'instances de la séquence représentée par ce seau.
lorsque tous les seaux ont été parcourus, 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 de vos limites. Si vous lisez les séquences à partir d'un fichier, au lieu de les conserver en mémoire, plus de mémoire serait disponible pour les seaux.

Je pense que cela serait plus rapide que de faire le tri in situ car les seaux ont probablement la bonne taille pour tenir dans votre ensemble de travail.

Voici une astuce qui montre la technique

#include 
#include 

#include 

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)
{
    // Trier cette séquence
    const char *testSequence = "CAGCCCAAAGGGTTTAGACTTGGTGCGCAGCAGTTAAGATTGTTT";

    size_t testSequenceLength = strlen(testSequence);

    setup();

    // charger les séquences dans les seaux
    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]++;
    }

    /* afficher les séquences triées */
    for (z = 0; z < bucketCount; z++)
    {
        while (bucket[z] > 0)
        {
            show(z);
            bucket[z]--;
        }
    }

    teardown();

    return 0;
}

Answer 5

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 éléments, int préfixe)
    if (éléments.Count < SEUIL)
         return InMemoryRadixSort(éléments, préfixe)
    else
         return DiskBackedRadixSort(éléments, préfixe)

DiskBackedRadixSort(éléments, préfixe)
    DiskBackedBuffer[] seaux
    foreach (élément in éléments)
        seaux[élément.MSB(préfixe)].Add(élément);

    List ret
    foreach (seau in seaux)
        ret.Add(sort(seau, préfixe + 1))

    return ret

Je recommande également d'essayer de regrouper dans un plus grand nombre de seaux, par exemple, si votre chaîne était :

GATTACA

le premier appel de MSB retournerait le seau pour GATT (256 seaux 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, donc expérimentez avec.