Transformation de Burrows-Wheeler


22 juillet 2007

J'avais besoin d'un programme qui réalisait la transformation de Burrows-Wheeler (TBW) sur de simples lignes de texte et non pas sur des énormes fichiers binaires comme c'est le cas d'habitude. J'ai donc pondu le code qui suit :

#include <iostream>

int LineLen;

int Cmp(const void* a,const void*b)
{
    return strncmp(*(const char**)a,*(const char**)b,LineLen);
}

int main()
{
    int i;
    char    *Line=new char[1000], **Rows=new char*[1000];
    while(gets(Line)){
        LineLen = strlen(Line);
        //On utilise un double buffer pour la ligne courante
        memcpy(Line+LineLen, Line, LineLen);
        //Les permutations en utilisant de simples pointeurs
        for(i=0; i<LineLen; i++)        Rows[i]=Line+i;
        //On trie les pointeurs en fonction de leurs données
        qsort(Rows, LineLen, 4, Cmp);
        //On cherche la ligne originale dans la liste triée.
        i=0;
        while (i<LineLen && strncmp(Rows[i], Line, LineLen))    i++;
        int OriginalIndex = i;
        //
        printf("%d\t",OriginalIndex);
        for(i=0; i<LineLen; i++)    printf("%c", Rows[i][LineLen-1]);
        printf("\n");
    }
    delete []Line;
    delete []Rows;
}

pour les données d'entrée suivante :

ABCABC
abracadabra
BurrowsWheelerTransform
gaattcctgggcctggggctgtggcagctgcctcgtccct
tcacctcctggcttattctctccctccatatcttagcaat
ctcatgcctgtaatcccagcattttggtaggccaaggcgg
gtcggatcacctgaggtcaggagttcgagggccagcctga
tgaccatggtgaaaccccatctctactaaaaatacaaaat
taatcgggcatggtggcacatgcctgtaatcccagctact
ctgaggcaggagaatcgcttaaacccaggaagtggaggtt
tcagtgagctgagattgtgccattgcactccagcctgggt
aacaagagcaaaactccatcaaaaaaaataatatatgtat
atatatattacaattttatatatatatacacattatgtaa
taccattttatatatatacattacgtaatggtaaatgttt
gatcgtctccctggagaataatccccaatgtgaaattact
ctaagtggtgggattacaggcgtgtgcccaacttttcctg
agccttttgaggctgacaccagaggtagaagcccagcctc
tccccactggccatgtggggagaggctccagcctgcagca
accagggatctggcctcaagtgatgccccaacagtgggcg
1001110010001111110101011101101110010100
1100111011000100110110010000011011100111
1000100001000001010101101000010110000011
0000000000111111000100100001011010100011
0011011000001111001010010001001011101111
0010110110101011100100010011000000010000
1101011110010100001000111001000100011001
0011001010001001001110100000101001001101
1110110100001001100110001010101001100100
1110001011101110110000110100011011100001

le programme renvoie :

0   CCAABB
2   rdarcaaaabb
0   mrsrhelsWerafreToruwnBo
16  gcagtgctgtccgccgtgtgagtggtgtctgtcccgccca
28  cctctatgtcttcatctctctgagttattcccctacccaa
17  ctctaacccctggctgggcagtggaacttgggcaatctta
31  cccgggggtctgagttcccttgggaacgaagaaagtgccg
37  tacgaaaaatgataccaacccacaatttttgcaccagtaa
31  ttctcaaccgcagctgtgcaggtatgcttggtcaaagacc
17  taggagccggaacgcattggagtggcaataaagtatcgcg
29  cgcgccggctctggagcttattaatgatgtccgtgctgaa
8   caaacaaatacataaagtaatcttgactaaataagaaaca
11  tacttctttttattttcaaaatgtaaaaataataaatata
25  ttatttttttaaccacaatgctggattaaataatatttga
25  ggtcatgaagaacccttcttactagttctatacaggccaa
13  cttacagcacgtggcagtgtgatttgaccttgcgggattc
7   ggctcaccggacctacgggcgcattaaaagaaggcctttc
35  cgccgcccgcgcctgctggcagatagaggtagttcaccga
3   ccagcacggctcacacgggctcgttgggtgataacagacg
21  1111100100111011001010010110111100010110
30  1010101011001100100011011111001110010010
30  1110100100000000001010111001001100000100
0   1000001010101000001010010101000100101110
8   1010101110010001010110010101010110111000
13  1010000010010010001010111100101010001000
35  1101110010011000010001000011101011010100
9   1010011101111110100000000000100010010100
39  1010101111110011000010011100010001000010
36  1110100000101001110111011110011010101000

Je voulais un code qui soit rapide à implémenter mais :

  • la limite de 1000 octets n'est pas terrible pour une routine généraliste.
  • j'utilise qsort qui n'est pas optimal par rapport à mon radix sort pour des données de taille importante.
  • J'utilise gets ce qui est déconseillé !
  • On n'a pas à dupliquer la ligne de données, en compliquant un peu le code ça devrait pouvoir marcher.
  • Le buffer Rows est de trop il faut trouver une astuce pour ne pas l'utiliser il prend une mémoire folle.

Tout ceci sera fixé plus loin, pour l'instant ce qui nous intéresse c'est d'avoir un outil de test à utiliser en conjonction avec d'autres compresseurs maison.

La transformation inverse

Maintenant que l'exécutable précédent (que j'ai appelé bwt) est fait et fonctionne, comment retrouver à partir de sa sortie les données originales ?

Le LZ* le plus intéressant : le LZSS


19 juillet 2007

J'apprécie depuis longtemps l'algorithme de compression sans perte LZSS et j'ai décidé de lui dédier une page de ce site pour l'illustrer et l'implémenter du mieux possible.

Pour faire vite et simple, l'algorithme compresse les données d'entrée en remplaçant des segments de caractères par des références vers des segments identiques déjà rencontrés. Ainsi Fulbert à la position 100 sera remplacé par Ful<25,4> si bert a déjà été rencontré à la position 25. Nous pouvons encoder la référence <position, taille> sur deux caractères, on verra plus tard comment et ainsi éliminer 4 - 2 = 2 caractères, d'où la compression. C'est, on s'en doute quelque chose qui fonctionne très bien pour le texte.

Voici venu le moment d'un exemple de code qui réalise cette compression dans des conditions idéales : on alloue un buffer suffisamment grand pour ne pas se poser la question du bouclage et on ne doit ni sauver les données ni les relire. Le but étant de tester le concept (qui marche depuis 1982 on ne fait donc pas là une grande découverte) et de fournir quelques statistiques :

getMatch(cursor)
    ... recherche ce qui se trouve à partir de la position <cursor> dans le dictionnaire.

store(position, len)
    ...stocke dans le dictionnaire <len> octets des données d'entrées à partir de <position>.

emitSegment(match)      output <match.offset, match.len>
emitChar(c)             output c

main()
{
    cursor = 0
    while (cursor != EOF)
        if (Match = getMatch(cursor))
            emitSegment(Match)
            store(cursor, Match.len)
            cursor += segment.length
        else
            emitChar(cursor)
            store(cursor, 1)
}

Le principe : on lit chaque caractère l'un après l'autre et soit on a déjà rencontré ce qui suit auquel cas on émet un segment, soit on n'a encore jamais rencontré ça et on émet le caractère tel quel. C'est extrèment simple et le code l'est également.

Tout ceci donne, une fois traduit, en C :

char*   Dic = NULL;
int     DicCursor = 0, DicLength = 0, DicStoredLength = 0;

void    addToDictionary(char c)
{
    Dic[DicCursor++] = c;
    DicStoredLength++;
}

bool    searchMatch(const char* strToCode, int cursor, int bufLn, int& longestMatchOffset, int& longestMatchLength)
{
    longestMatchLength = longestMatchOffset = 0;
    int DicPos = 0, i = 0;
    while (DicPos < DicStoredLength)        {
        while (DicPos+i<DicStoredLength && cursor+i<bufLn && strToCode[cursor+i] == Dic[DicPos+i])      i++;
        if (longestMatchLength<i)       {
            longestMatchOffset = DicPos;
            longestMatchLength = i;
        }
        DicPos++;
        i = 0;
    }
    return longestMatchLength != 0;
}


void    outputChar(char c)              {   printf("%c", c);    }
void    outputCode(int offset, int ln)  {   printf("<%d, %d>", offset, ln);     }

int main(int argc, char *argv[])
{
    DicLength = 2048;
    FILE*   FileToCode = fopen("../../config.h.in", "r");
    if (!FileToCode)        return EXIT_FAILURE;
    fseek(FileToCode, 0, SEEK_END);
    long    FileSize = ftell(FileToCode);
    char*   ToCode = new char[FileSize];
    rewind(FileToCode);
    fread(ToCode, 1, FileSize, FileToCode);
    fclose(FileToCode);
    //
    Dic = new char[DicLength];
    int Ln = FileSize;
    int i = 0, EncodedSize = 0;
    int NbRawChars = 0, NbBlocks = 0;
    int HeaderSize = 3, j;
    while (i<Ln)        {
        int MatchOffset, MatchLength;
        char    c = ToCode[i];
        if (searchMatch(ToCode, i, Ln, MatchOffset, MatchLength) && MatchLength>2)      {
            outputCode(MatchOffset, MatchLength);
            for(j=i; j<i+MatchLength; j++)      addToDictionary(ToCode[j]);
            i += MatchLength;
            EncodedSize += 2;
            NbBlocks++;
        }   else    {
            outputChar(c);
            i++;
            EncodedSize++;
            NbRawChars++;
            addToDictionary(c);
        }
    }
    delete []Dic;
    printf("\n");
    printf("Taille originale : %d\n", Ln);
    printf("Taille compressée hors header et flags : %d\n", EncodedSize);
    printf("Nombre de caractères émis de façon brute : %d\n", NbRawChars);
    printf("Nombre de blocs <offset,ln> : %d\n", NbBlocks);
    int TotalPackedSize = EncodedSize + HeaderSize + (NbRawChars>>3) + 1;
    printf("Taille totale compressée : %d\n", TotalPackedSize);
    printf("Ratio : %.3f %%\n", 100.0f * TotalPackedSize / Ln);
    return EXIT_SUCCESS;
}

Voyons ce que ça donne pour plusieurs sources d'entrée :

Taille originale : 35147
Taille compressée hors header et flags : 11621
Nombre de caractères émis de façon brute : 1933
Nombre de blocs <offset,ln> : 4844
Taille totale compressée : 11866
Ratio : 33.761 %", "Source : <a href='http://www.gnu.org/licenses/gpl-3.0.txt'>GPL v3</a>
Taille originale : 702748
Taille compressée hors header et flags : 174968
Nombre de caractères émis de façon brute : 3628
Nombre de blocs <offset,ln> : 85670
Taille totale compressée : 175425
Ratio : 24.963 %

Vous remarquerez que la taille compressée est inférieure de 80ko au fichier Zip référencé sur le site du Projet Gutenberg. C'est un bon début :)

Bien sûr tout ça n'est pas très sérieux : ça prend un temps fou (la méthode de comparaison est risible, au moins 3 minutes pour Pride and Prejudice...), le buffer s'adapte aux données d'entrées (je préfèrerai qu'il soit fixe) et il n'y a toujours pas de sauvegarde des données pour vérifier que ça fonctionne dans les deux sens. Mais ça n'est pas grave, c'est la première version et on a vérifié que le code est faisable (ouf c'était dur) (je plaisantais dans la dernière parenthèse, on est d'accord hein ?).

La première optimisation

En rajoutant trois petites lignes à la recherche on divise par 3 le temps passé dans la routine. Ces lignes sont là pour passer sur un caractère s'il n'est pas égal au premier recherché :

bool    searchMatch(const char* strToCode, int cursor, int bufLn, int& longestMatchOffset, int& longestMatchLength)
{
    longestMatchLength = longestMatchOffset = 0;
    //Recherche naïve
    int DicPos = 0, i = 0;
    while (DicPos < DicStoredLength)        {
        char    CharToLookFor = strToCode[cursor];
        while (DicPos<DicStoredLength &&  Dic[DicPos] != CharToLookFor)     DicPos++;
        i=1;
        while (DicPos+i<DicStoredLength && cursor+i<bufLn && strToCode[cursor+i] == Dic[DicPos+i])      i++;
        if (longestMatchLength<i)       {
            longestMatchOffset = DicPos;
            longestMatchLength = i;
        }
        DicPos++;
        i = 0;
    }
    return longestMatchLength != 0;
}

Deuxième optimisation : hash table et listes chaînées

Tout le monde connait le concept et l'intérêt des Hash Table, et bien pour le LZSS c'est impressionnant d'intérêt. Si on cherche une chaîne dans le dictionnaire on peut parcourir ce dernier comme précédemment ou bien demander à une fonction "donne-moi les emplacements du dictionnaire qui contiennent telle lettre". Plus de parcours ! On va directement au bon endroit. La lettre en question s'appelle la clé (de la requête).

Bien. Mais quelle est la structure de donnée qui permet de poser cette question efficacement, comprendre rapidement (en O(1) si possible) ? Et bien c'est une Hash Table, dont les clés sont les premières lettres hashées, combinée avec des listes chaînées contenant des pointeurs vers les emplacement dans le dictionnaire.

Tout ça économise énormément de temps sur la recherche d'un segment du dictionnaire intéressant mais c'est optimisable : il suffit non plus d'utiliser la première lettre mais plusieurs :) En fait je suis arrivé à une version qui utilise les 4 premières lettres des données d'entrées, les combine en un nombre qui représente un index dans ma Hash Table pour parcourir une liste de cellules chaînées pour trouver mon segment le plus intéressant. Les résultats ?

Un résultat intermédiaire m'a donné pour Pride & Prejudice 102 secondes pour une recherche naïve (qui utilise tout de même la première optimisation) et 3.37sec pour la version "Hash Table" avec un buffer illimité.

Si on est maintenant prêt à oublier un peu le taux de compression pour privilégier la vitesse d'encodage, en limitant la taille du dictionnaire à 32768 octets on passe à 13.7 secs pour la version naïve et à 0.27 secs pour la version utilisant une hash table.

Rien à dire ça commence à être praticable...

Reste à présent à implémenter la sérialisation (la sauvegarde et la relecture) et à améliorer l'occupation mémoire parce que pour arriver à presque 1/4 de seconde j'ai créé plein de choses dans la RAM du PC...

La moralité : quand on constate que j'ai un facteur 1000 entre cette optimisation et la recherche naïve on se dit que l'algo le plus intéressant de recherche de chaîne de caratères ne pourra jamais battre une indexation à 4 caractères ! Dommage, parce que ça veut dire que Jalut écrase Dāwūd, mais c'est logique. D'un autre côté si quelqu'un est arrivé à faire quelque chose avec KMP ou Boyer-Moore ça m'intéresse !

La suite...

on décode des données sauvées On compare toujours les temps des différentes méthodes pour savoir où l'on en est. L'usage de l'algo ? On utilise des fichiers de différents types : wav, dwg, txt & cpp, images raw. Les améliorations : - taille variable du stockage des segments (1, 2 ou 3 octets). - on sauve séparément les flux de caractères, de bits de contrôle et de segments (position et taille séparément également). - on compresse à nouveau les flux différemment en fonction de leurs contenus - stockage définitifs des segments à partir d'un certain seuil (ils ne sont plus supprimés du dico mais déplacés ailleurs)

Accueil