Problème de consommation mémoire

Version imprimable

Bonjour à tous,

Voilà je travail actuellement sur un programme qui analyse un graphe de séquences d'ADN en entrée et effectue des calculs dessus.
Comme l'ADN ne comporte que 4 nucléotides ('A', 'C', 'G' et 'T') j'ai développé une suite de méthodes permettant de stocker 3 nucléotides sur un octet (type unsigned char) afin de gagner de la mémoire. Une séquence est donc représenté par un tableau de unsigned char.
Le système marche comme cela : chaque nucléotide utilise 2 bits de l'octet ( A = 00, C = 01, G = 10, T=11). Les 6 bits de poids fort (les plus a gauche) de l'octet codent donc 3 nucléotides. Les 2 bits restants indique le décalage (dans le cas ou un char ne stock pas 3 caractères) : 11 = ce char contient 3 nucléotides, 10 = ce char ne contient que deux nucléotides, 01 = ce char ne contient qu'une nucléotide. Pour des raisons lié aux structures que j'utilise dans mon programme, l'option "ce char contient 3 nucléotides" n'est pas représenté par 0 car la séquence "AAA" sera alors représenté par 00000000, qui est le caractère \0. Enfin, toute séquence se termine par le caractère '\0' comme dans un vrai string.

J'ai écrit deux versions de mon programme, l'une avec ce système et l'autre sans ce système (où donc une nucléotide = un char et une séquence = un tableau de char). Mon problème est le suivant : le programme qui utilise de simples string consomme moins de mémoire que celui qui utilise mon système pour coder 3 nucléotides sur un octet :mrgreen:

Je tiens a préciser quelques petites choses :
- Mes deux programmes sont certifiés 0 fuites mémoires par Valgrind
- Pour mesurer la consommation mémoire, j'ai utilisé la commande smem au sein de mon programme et l'outil tstime, en regardant à chaque fois le RSS (Resident Set Size).
- Je ne suis pas un grand spécialiste du C :D

Voici le code de mon système (system.c) :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "sequence.h"
#include <math.h>
 
#define MASK_0 3
#define MASK_1 12
#define MASK_2 48
#define MASK_3 192
#define MASK_1_2_3 252
#define MASK_0_1_2_3 255
 
unsigned char* createSequence(char *str){
    int lg, dec, lg_str, current_dec;
    char* seq;
 
        lg_str = strlen(str);
        dec = lg_str%3;
 
        //Calcul de la longueur du tableau de caractère
        if (dec == 0) lg = lg_str/3;
        else lg = lg_str/3+1;
 
        //Initialisation du tableau de caractère
        seq = calloc(lg+1, sizeof(unsigned char));
 
        if (seq != NULL){
 
            if (dec == 0) dec = 3;
            current_dec = 3;
 
            //On parcours le chaine en entrée et tous les 4 caractères, on insère dans le tableau après tranformation.
            int i, k;
            for (i=0; i<lg; i++){
                if (i==lg-2) current_dec = dec;
 
                for (k=0; k<current_dec; k++){
                    switch(str[i*3+k]){
                        case 'a':
                        case 'A': seq[i] = (seq[i]<<2) | 0;
                                  break;
                        case 'c':
                        case 'C': seq[i] = (seq[i]<<2) | 1;
                                  break;
                        case 'g':
                        case 'G': seq[i] = (seq[i]<<2) | 2;
                                  break;
                        case 't':
                        case 'T': seq[i] = (seq[i]<<2) | 3;
                                  break;
                        case ' ':
                        case '\n':
                        case '\0': break;
                        default : printf("Invalid character encountered : - Operation aborted");
                                  exit(0);
                    }
 
                    if (current_dec != 3) seq[i] = (seq[i]<< (8-dec*2));
                }
 
                seq[i] = (seq[i]<<2) | MASK_0;
            }
            seq[i-1] = ((seq[i-1] & MASK_1_2_3) << (6-dec*2)) | dec;
            seq[lg] = '\0';
 
            return seq;
        }
        else {
            fprintf (stderr, "Memoire insufisante\n");
            exit (EXIT_FAILURE);
        }
}
 
int seqlen(unsigned char* seq){
    int lg = strlen(seq);
    int dec = seq[lg-1] & MASK_0;
 
    if (dec == 3) return lg*3;
    return (lg-1)*3+dec;
}
 
int seqcmp(unsigned char* seq1, int position1, unsigned char* seq2, int position2, int costMismatch){
 
    int offset2=position2%3, offset1 = position1%3;
    int nb1, nb2;
 
 
    switch (offset1){
        case 0: nb1 = seq1[position1/3] & MASK_3;
        case 1: nb1 = seq1[position1/3] & MASK_2;
        case 2: nb1 = seq1[position1/3] & MASK_1;
        default: return -1;
    }
 
    switch (offset2){
        case 0: nb2 = seq2[position2/3] & MASK_3;
        case 1: nb2 = seq2[position2/3] & MASK_2;
        case 2: nb2 = seq2[position2/3] & MASK_1;
        default: return -1;
    }
 
    if (offset2>offset1){
        if (nb2 == pow(3,offset2-offset1)*nb1) return 0;
        return costMismatch;
    }
    if (nb2 == pow(3,offset1-offset2)*nb1) return 0;
    return costMismatch;
}
 
unsigned char* reverseSequence(unsigned char* seq, int complement){
    int i, lg, dec;
    int seq_lg = seqlen(seq);
 
    //Calcul de la longueur du tableau de caractère
    if (seq_lg%3 == 0) lg = (seq_lg/3);
    else lg = (seq_lg/3)+1;
 
    dec = seq_lg%3;
    if (dec == 0) dec = 3;
 
    char *newSeq = calloc(lg+1, sizeof(unsigned char));
 
    if (newSeq != NULL){
        if (complement==0){
            for (i=0; i<lg; i++){
                newSeq[i] = ~(seq[i]) & MASK_1_2_3 | MASK_0;
            }
            newSeq[i-1] = (newSeq[i-1] & (MASK_0_1_2_3 << 8-dec*2)) | dec;
            newSeq[lg] = '\0';
            return newSeq;
        }
 
        for (i=0; i<lg/2; i++){
            newSeq[i] = reverseCompChar(seq[lg-i-1]);
            newSeq[lg-i-1] = reverseCompChar(seq[i]);
        }
 
        if (lg%2 == 1) newSeq[i] = reverseCompChar(seq[i]);
 
        newSeq[lg-1] = newSeq[lg-1] | MASK_0;
        newSeq[lg] = '\0';
 
        if (seq_lg%3 != 0){
            char* res = subSequence(newSeq, 3-(seq_lg%3), lg*3-1);
            free(newSeq);
            return res;
        }
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* subSequence(unsigned char* seq, int posDeb, int posFin){
 
    int i, decTmp = (posDeb%3)*2, k=0, lg, newSeq_lg, seq_lg, dec;
    unsigned char tmp1;
    char* newSeq;
 
    //Calcul de la longueur du tableau de caractère
    newSeq_lg = posFin-posDeb+1;
    seq_lg = seqlen(seq);
    dec = newSeq_lg%3;
 
    if (dec == 0) lg = newSeq_lg/3;
    else lg = newSeq_lg/3+1;
    newSeq = calloc(lg+1,sizeof(unsigned char));
 
    if (newSeq != NULL){
        for (i = posDeb/3; i<posFin/3; i++){
            tmp1 = (seq[i] & MASK_1_2_3) << decTmp;
            if (seq_lg > 3) tmp1 = tmp1 | (seq[i+1] >> (6-decTmp));
            newSeq[k] = tmp1 & MASK_1_2_3 | MASK_0;
            k++;
        }
 
        if (posDeb%3!=0) newSeq[k] = (seq[i] & MASK_1_2_3) << decTmp;
        else{
            if ((lg != 1) || (lg==1 && ((posDeb/3)==(posFin/3)))){
                if (dec == 0) newSeq[k] = (seq[i] & MASK_1_2_3)/* & MASK_0_1_2_3*/;
                else newSeq[k] = (seq[i] & MASK_1_2_3) & (MASK_0_1_2_3 << 8-dec*2);
            }
        }
 
        if (dec == 0) dec = 3;
        if (dec!=0) newSeq[lg-1] = newSeq[lg-1] | dec;
        newSeq[lg] = '\0';
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* concatSequence(unsigned char* seq1, unsigned char* seq2){
    int lgTab1, lgTab2, i, k, decTmp, lgTmp;
    int seq_lg, seq1_lg, seq2_lg;
    char* seq;
 
    // Oninitialise la nouvelle séquence
    seq1_lg = seqlen(seq1);
    seq2_lg = seqlen(seq2);
 
    seq_lg = seq1_lg + seq2_lg;
 
    if (seq_lg%3 == 0) lgTmp = seq_lg/3;
    else lgTmp = (seq_lg/3)+1;
 
    // Calcul de la longueur du tableau de la séquence 1
    if (seq1_lg%3 == 0) lgTab1 = seq1_lg/3;
    else lgTab1 = (seq1_lg/3)+1;
 
    // Calcul de la longueur du tableau de la séquence 2
    if (seq2_lg%3 == 0) lgTab2 = seq2_lg/3;
    else lgTab2 = (seq2_lg/3)+1;
 
    seq = calloc(lgTmp+1, sizeof(unsigned char));
 
    // Recopie du tableau de la séquence 1 dans la nouvelle séquence
    for (i=0; i<lgTab1; i++){
        seq[i] = seq1[i];
    }
 
    decTmp = (seq1[i-1] & MASK_0) * 2;
    seq[i-1] = seq[i-1] & MASK_1_2_3 | MASK_0;
 
    // Si aucun décalage dans la séquence 1, je recopie le tableau de la séquence 2 à la suite du 1 dans la nouvelle séquence
    if (decTmp == 6){
        for (k=0; k<lgTab2; k++){
            seq[i] = seq2[k];
            i++;
        }
    }
    else {
        // Gestion du décalage entre la dernière ligne du premier tableau et la première ligne du dernier tableau
 
        seq[i-1] = seq[i-1] | (seq2[0] >> decTmp) | MASK_0;
        if (i < lgTmp) seq[i] = (seq2[0] & MASK_1_2_3) << (6-decTmp);
 
        // Recopie du deuxième tableau avec gestion du décalage
        for (k=0; k<lgTab2-1; k++){
            seq[i] = seq[i] | (seq2[k+1] >> decTmp) | MASK_0;
            seq[i+1] = (seq2[k+1] & MASK_1_2_3) << (6-decTmp);
            i++;
        }
 
        if (seq_lg%3 == 0) seq[lgTmp-1] = seq[lgTmp-1] | MASK_0;
        else seq[lgTmp-1] = seq[lgTmp-1] | (seq_lg%3);
    }
 
    seq[lgTmp]= '\0';
 
    return seq;
}
 
char *substr(const char *src,int pos,int len) {
  char *dest=NULL;
  if (len>0) {
    dest = calloc(len+1, sizeof(char));
    if(NULL != dest) {
        strncat(dest,src+pos,len);
    }
  }
  return dest;
}
 
unsigned char reverseCompChar(unsigned char c){
    int k;
    unsigned char tmp = 0;
 
    for (k=0; k<3; k++){
        switch(k){
            case 0: tmp = ((~(c >> 2) & 3) | tmp) << 2; break;
            case 1: tmp = ((~(c >> 4) & 3) | tmp) << 2; break;
            case 2: tmp = ((~(c >> 6) & 3) | tmp) << 2; break;
        }
    }
 
    tmp = tmp | MASK_0;
 
    return tmp;
}

Puis le header (sequence.h) :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
#ifndef DEF_SEQUENCE
#define DEF_SEQUENCE
 
unsigned char* createSequence(char *str);
int seqlen(unsigned char* seq);
int seqcmp(unsigned char* seq1, int position1, unsigned char* seq2, int position2, int costMismatch);
unsigned char* reverseSequence(unsigned char* seq, int complement);
unsigned char* subSequence(unsigned char* seq, int posDeb, int posFin);
unsigned char* concatSequence(unsigned char* seq1, unsigned char* seq2);
char *substr(const char *src,int pos,int len);
unsigned char reverseCompChar(unsigned char c);
 
#endif

J'ai mesuré la mémoire à différents points de mon programme et rien que le début, n'utilisant que createSequence() utilise déjà beaucoup trop de mémoire. Voici un exemple d'appel a createSequence():

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
 
void debut_element(void *user_data, const xmlChar *name, const xmlChar **attrs) {
    DataParser* dp = (DataParser*)user_data;
     if (xmlStrEqual(name, BAD_CAST "node")) {
         Node* n = NULL;
         char* id = NULL;
         char* att = NULL;
         if (NULL != attrs) {
            int i;
            for (i = 0; attrs[i] != NULL; i+=2) {
                if (xmlStrEqual(attrs[i], BAD_CAST dp->nodeAttribute)) {
                    att = calloc(strlen((char *) attrs[i+1]) + 1, sizeof(char));
                    strcpy(att, (char *) attrs[i+1]);
                }
                else if (xmlStrEqual(attrs[i], BAD_CAST "id")) {
                    id = calloc(strlen((char *) attrs[i+1]) + 1, sizeof(char));
                    strcpy(id, (char *) attrs[i+1]);
                }
            }
 
            if (id != NULL){
                if (att != NULL){
                    dp->nb_nucl += (long) strlen((char *) att);
 
                    n = createNode(id, att);
                    int ret;
                    khiter_t k;
 
                    k = kh_put(32, dp->net, id, &ret);
                    kh_value(dp->net, k) = n;
 
                    free(att);
                }
            }
            else {
                fprintf (stderr, "Missing element(s) on one node while parsing graph\n");
                exit (EXIT_FAILURE);
            }
         }
     }
     else if (xmlStrEqual(name, BAD_CAST "edge")) {
        if (NULL != attrs) {
            int i;
            char* reading = NULL;
            khiter_t s, t;
            for (i = 0; attrs[i] != NULL; i+=2) {
                if (xmlStrEqual(attrs[i], BAD_CAST dp->edgeAttribute)) {
                    int l = strlen((char *) attrs[i+1]) + 1;
                    reading = calloc(l, sizeof(char));
                    strcpy(reading, (char *) attrs[i+1]);
                }
                else if (xmlStrEqual(attrs[i], BAD_CAST "source")) {
                   s = kh_get(32, dp->net, attrs[i+1]);
                }
                else if (xmlStrEqual(attrs[i], BAD_CAST "target")) {
                   t = kh_get(32, dp->net, attrs[i+1]);
                }
            }
 
            if (reading != NULL && s && t){
                addNeigh((Node*)kh_value(dp->net,s), (Node*)kh_value(dp->net,t), reading);
            }
            free(reading);
         }
     }
}

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
Node* createNode(char* id, char* attribute){
    Node* n = malloc(sizeof(Node));
    if ( n!= NULL ){
        n->id = id;
        n->attribute = createSequence(attribute);
        n->children = list_create();
        n->relative_position = list_create();
        return n;
    }
    else {
        fprintf (stderr, "insufficient memory to create a node\n");
        exit (EXIT_FAILURE);
    }
}

Voilà si quelqu'un à ne serait-ce qu'une idée, une piste, quelque chose, je lui en serait super reconnaissant ! Merci ;)

-createSequence() ne marche pas correctement si on n'a pas un nombre multiple de 3 de ACGT (des confusions dans le code entre dec et current_dec entre autres erreurs).

- De plus, il ne fonctionnera pas si il y a des blancs ou des '\n' dans le fichier. Le résultat sera faux (on traite seq[i] même dans ce cas : if (current_dec != 3) seq[i] = (seq[i]<< (8-dec*2)); .

- De toute façon, createSequence() alloue trop de mémoire dans ce cas, puisque strlen() compte les blancs et les '\n'.

- Le code ci-dessous essaye de pallier à ces défauts (printACGT() utilitaire pour visualiser ce qu'on obtient):

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
enum {NUCLEO_A, NUCLEO_C, NUCLEO_G, NUCLEO_T, NUCLEO_IGN};
//--------------------------------------------------------------
unsigned char* createSequence(char *str)
{
    int lg, i, dec;
    int nucleo;
    char * seq;
    char * pseq;
    //Initialisation du tableau de caractère
    lg = strlen(str);
    dec = lg%3;
    if (dec == 0) lg = lg/3;
    else lg = lg/3+1;
    seq = calloc(lg+1, sizeof(unsigned char));
    if (seq != NULL)
    {
      //lire la séquence
       for(i = 0, dec = 0 ;*str != '\0'; str++)
       {
          switch(*str)
          {
             case 'a' :
             case 'A' : nucleo = NUCLEO_A;   break;
             case 'c' :
             case 'C' : nucleo = NUCLEO_C;   break;
             case 'g' :
             case 'G' : nucleo = NUCLEO_G;   break;
             case 't' :
             case 'T' : nucleo = NUCLEO_T;   break;
             case ' ' :
             case '\n': nucleo = NUCLEO_IGN; break;
             default  : printf("Invalid character encountered : - Operation aborted");
                        free(seq);
                        exit(0);
          }
          if(nucleo != NUCLEO_IGN) // ne traiter que si ACGT
          {
             seq[i] = (seq[i]<<2) | nucleo;
             dec++;
             if(dec==3) // octet complet
             {
                seq[i] = (seq[i]<<2)+3;
                i++;
                dec= 0;
             }
          }
       } // fin de lecture
       if(dec != 0)  // si le dernier octet est incomplet
       {
          seq[i] = (seq[i]<< (8-dec*2))+dec;
          i++;
       }
       pseq = realloc(seq,i+1);    // redimensionner le tableau
       if(pseq != NULL) seq = pseq;
       seq[i] = '\0';             // et le terminer
    }
    return seq;
}
//--------------------------------------------------------------
void printACGT( unsigned char c)
{
  static char t[] = "ACGT";
  int count = c&3;
  printf(" nb = %d (%02x) -> ",count,c);
  if(count >=1) printf("%c", t[(c&0xc0)>>6]);
  if(count >=2) printf("%c", t[(c&0x30)>>4]);
  if(count >=3) printf("%c", t[(c&0x0C)>>2]);
  printf("\n");
}
//--------------------------------------------------------------
int main(void)
{
     unsigned char * p =  createSequence("  a c gtcc ggcc \n ggtt ");
     unsigned char *pp;
     if(p != NULL)
     {
        for(pp=p; *pp != 0; pp++)printACGT(*pp);
        free(p);
     }
     else
     {
       fprintf (stderr, "Memoire insufisante\n");
       return EXIT_FAILURE;
     }
     return EXIT_SUCCESS;
}
//---------------------------------------------------------------------------

30/07/2012, 09h14
boboss123

salut,

J'ai quelques remarques :
- Pourquoi mettre un '\0' en fin d'un tableau qui contient des données "non ascii" (binaire) ?
- Pourquoi veux-tu coder tes séquences sur des unsigned char ? si c'est pour accélérer le rapidité des tes calculs, je pense qu'il y a 95% de chance que ça va plutôt te les ralentir. L'avantage que je vois pour cette méthode, c'est si tu veux stocker tes données sur un disque ou que tu n'as pas assez de RAM.
30/07/2012, 11h17
Obsidian

Citation:

Envoyé par boboss123

salut,

J'ai quelques remarques :
- Pourquoi mettre un '\0' en fin d'un tableau qui contient des données "non ascii" (binaire) ?

Ben pour connaître la fin du tableau. Comme pour une chaîne.

Citation:

- Pourquoi veux-tu coder tes séquences sur des unsigned char ? si c'est pour accélérer le rapidité des tes calculs, je pense qu'il y a 95% de chance que ça va plutôt te les ralentir.

Pourquoi ?

Citation:

L'avantage que je vois pour cette méthode, c'est si tu veux stocker tes données sur un disque ou que tu n'as pas assez de RAM.

Il me semble que ce qui le préoccupe dès le départ, c'est ça consommation mémoire. Il faut savoir qu'en bioinformatique, quand on travaille sur des séquences d'ADN, on manipule un très grand nombre de très grandes chaînes. Le chromosome 1 humain à lui seul représente 250 millions de paires de bases. Tu le découpes ensuite en n fragments que tu répliques au moins une trentaine de fois pour faire des statistiques, et tu comprends pourquoi on essaie d'être économe.

Citation:

Envoyé par Obsidian

Pourquoi ?

Par ce que s'il faut faire des calculs (ex : compter le nombre de Nucleo A), il faudra faire des masques avant de faire la comparaison. Un tableau de int sera plus optimale pour faire des calculs il me semble (-1 = fin du tableau).
=> mais vu que je ne sais pas exactement quelle sorte de calculs il fait, c'est à confirmer...
=> s'il y a un manque de RAM, peut être que le traitement de la chaine peut être découpé
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 int main() { /* char tab[100000]; char * ptchar; char val; // */ //* // c'est plus rapide int tab[100000]; int * ptchar; int val; // */ int i,j; for(j=0; j<100; j++){ for(i=0; i<100; i++){ ptchar = tab; unsigned int len = sizeof(tab) / sizeof(*tab); while(len-- > 0){ val += *ptchar++; } } } return 0; }
Aussi pour la détection de la fin de la chaine un while(len--){} est peut être plus rapide qu'un while(val != 0x00){} => faut voir l'algo en entier.

=> je voulais juste mettre en évidence que ce n'est pas parce qu'une application utilise moins de RAM, que le traitement va aller plus vite

30/07/2012, 14h29
Obsidian
Citation:

Envoyé par boboss123

=> je voulais juste mettre en évidence que ce n'est pas parce qu'une application utilise moins de RAM, que le traitement va aller plus vite

Non, en effet. C'est même généralement le contraire en temps normal, mais ce n'est pas une relation linéaire, loin de là. C'est vraiment à examiner au cas par cas, mais c'est une pratique qui tend à se perdre, malheureusement.

Citation:

Par ce que s'il faut faire des calculs (ex : compter le nombre de Nucleo A), il faudra faire des masques avant de faire la comparaison. Un tableau de int sera plus optimale pour faire des calculs il me semble (-1 = fin du tableau). => mais vu que je ne sais pas exactement quelle sorte de calculs il fait, c'est à confirmer...

Ce qui est tueur en terme de performances, c'est le désalignement des données, spécialement lorsque tu effectues des accès aléatoires. Pour le reste, ça se justifie lorsque la quantité de données reste raisonnable et que le domaine des valeurs possibles occupe une place comparable à celle d'un int.

Par contre, dans le cas précis, il s'agit quand même d'utiliser 32 bits — voire 64 — pour stocker des éléments de 2 bits ! Soit 16 à 32 fois la mémoire réellement nécessaire. Et comme on manipule couramment, dans le domaine qui occupe présentement, des séquences de plusieurs centaines de méga-octets, le surcoût en mémoire cesse d'être négligeable.

En plus, en l'occurrence, les tailles des entiers usuels sont exactement multiples de la taille de l'élément à stocker, multiples qui sont eux-mêmes des puissances de 2. Si bien que les offsets peuvent se calculer avec de simples décalages de bits.

Il faut également tenir compte du fait que ces séquences ont aussi vocation à être enregistrés sous de nombreuses version sur le disque, ainsi qu'à être transmises par le réseau.

Enfin, le comptage des « A » ou autre valeur est le mauvais exemple puisque c'est typiquement une opération de lecture linéaire, ce qui est le meilleur cas en terme de complexité. Donc, avec des nucléotides empaquetés de la sorte, avec un seul accès bus, tu charges d'un coup 24 nucléotides en 32 bits et 48 en 64 bits. Et le comptage, non seulement ne fait plus appel au bus, mais n'utilise plus non plus la ligne de cache ! Tout reste en interne à l'intérieur du processeur jusqu'au bloc suivant.
Citation:
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 int main() { /* char tab[100000]; char * ptchar; char val; // */ //* // c'est plus rapide int tab[100000]; int * ptchar; int val; // */ …
Chez moi, c'est le contraire (Core 2 Quad Q9550 64 bits ; 4 Gio ; Fedora 16 ; GCC 4.6.3). Quand je fais des « time » successifs, j'obtiens les résultats suivants :
char :
Code:

1 2 3 4 5 6 7 real 0m2.495s real 0m3.287s real 0m2.499s real 0m2.489s real 0m3.209s real 0m3.026s real 0m2.497s
int :
Code:

1 2 3 4 5 6 7 real 0m3.354s real 0m3.186s real 0m2.493s real 0m3.305s real 0m3.324s real 0m3.277s real 0m3.205s
Citation:

Aussi pour la détection de la fin de la chaine un while(len--){} est peut être plus rapide qu'un while(val != 0x00){} => faut voir l'algo en entier.

Et comment fais-tu pour calculer « len » si tu ne disposes pas à l'avance de cette information ?

En ce qui concerne la différence entre « len-- » et « val != 0x00 », le compilateur va lui-même procéder à des optimisations qui rendront ces deux versions indistinguables dans le code final. Mais en admettant que ce ne soit pas le cas, à partir du moment où tu lis « val » à l'intérieur de la boucle, ça revient exactement au même :

« len-- » deviendra faux lorsqu'il arrivera à zéro, et le micro-processeur se sert d'un détecteur de zéro pour le savoir (un NON-OU électronique sur chaque bit du registre dont le résultat correspond au bit Z des flags). Or, une « comparaison » (avec CMP en assembleur) est faite par le micro-processeur en soustrayant les opérandes et l'égalité est établie … si le résultat est nul ! Donc en examinant le même bit Z.

Comme, en plus, tu cherches à savoir si « var » est ou non différente de zéro, il n'y a même pas d'opération à effectuer : le simple fait d'avoir chargé sa valeur en mémoire pour l'exploiter te donne la réponse dans Z.
30/07/2012, 14h50
boboss123

Citation:

Envoyé par Obsidian

Chez moi, c'est le contraire (Core 2 Quad Q9550 64 bits ; 4 Gio ; Fedora 16 ; GCC 4.6.3). Quand je fais des « time » successifs, j'obtiens les résultats suivants :

Tu as une explication pour ça (int fait bien 64bits ? tu as bien un OS 64 bits ?) ?

Ha messieurs, je suis content de savoir que ce sujet amène à discussion :)

Je voudrais d'abord remercier diogène pour son code : avec celui-ci maintenant, mes deux versions de mon programme utilisent dans tous les cas ... exactement la même quantité de mémoire :mrgreen: Mais je suppose qu'arriver à ce stade, il y a quelque chose qui cloche non pas dans le système que je vous présente mais ailleurs dans le programme ... Vous trouverez en bas de ce post la nouvelle version de sequence.c corrigé (il subsistait des erreurs dans subSequence() et concatSequence() mais rien qui influence la taille de la mémoire).

Pour répondre à boboss123 et Obsidian, effectivement je travail dans le domaine de la bioinformatique et la taille des données en entrée de ce programme peut se compter en Tera-octet. Il y a un choix a faire entre un programme optimal niveau temps et un programme optimal niveau mémoire et en bioinfo, le choix s'impose de lui même : si je n'ai pas assez de mémoire pour analyser mon fichier en entrée, il ne me sert à rien de gagner du temps.
Pour répondre à boboss123, effectivement un système à base de tableau d'int aurait été plus simple. Avec ce système, j'aurais pu coder 4 nucléotides sur un octet, et me servir d'un octet de fin pour indiquer à la fois la fin de la séquence et le décalage dans le dernier octet. J'avais même implémenté ce système :) Le problème, c'est que mes séquences d'ADN encodé de cette manière servent de clés dans une HashTable qui se sert des fonctions string pour comparer, copier, etc ... les clés. Pas de '\0' a la fin et toutes ces fonctions ont des comportements anormaux. Si j'utilise ce système avec des chars, il suffit d'un seul octet sensé représenter AAAA et j'obtiens le caractères '\0' signalant la fin de ma chaine, prématurément. C'est là que vous me direz : avec les fonctions mem, ça aurait marché si on connait la taille des clés à l'avance. Je l'ai cru avec la librairie Google SparseHash qui gère les fonctions mem si on lui précise mais malgré tout, j'obtenais des bugs énormes. J'ai donc changé de librairie et je me suis tourné vers Khash, très efficace.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "sequence.h"
#include <math.h>
 
#define MASK_0 3
#define MASK_1 12
#define MASK_2 48
#define MASK_3 192
#define MASK_1_2_3 252
#define MASK_0_1_2_3 255
 
enum {NUCLEO_A, NUCLEO_C, NUCLEO_G, NUCLEO_T, NUCLEO_IGN};
//--------------------------------------------------------------
unsigned char* createSequence(char *str)
{
    int lg, i, dec;
    int nucleo;
    unsigned char * seq;
    unsigned char * pseq;
    //Initialisation du tableau de caractère
    lg = strlen(str);
    dec = lg%3;
    if (dec == 0) lg = lg/3;
    else lg = lg/3+1;
    seq = calloc(lg+1, sizeof(unsigned char));
    if (seq != NULL)
    {
      //lire la séquence
       for(i = 0, dec = 0 ;*str != '\0'; str++)
       {
          switch(*str)
          {
             case 'a' :
             case 'A' : nucleo = NUCLEO_A;   break;
             case 'c' :
             case 'C' : nucleo = NUCLEO_C;   break;
             case 'g' :
             case 'G' : nucleo = NUCLEO_G;   break;
             case 't' :
             case 'T' : nucleo = NUCLEO_T;   break;
             case ' ' :
             case '\n': nucleo = NUCLEO_IGN; break;
             default  : printf("Invalid character encountered : - Operation aborted");
                        free(seq);
                        exit(0);
          }
          if(nucleo != NUCLEO_IGN) // ne traiter que si ACGT
          {
             seq[i] = (seq[i]<<2) | nucleo;
             dec++;
             if(dec==3) // octet complet
             {
                seq[i] = (seq[i]<<2)+3;
                i++;
                dec= 0;
             }
          }
       } // fin de lecture
       if(dec != 0)  // si le dernier octet est incomplet
       {
          seq[i] = (seq[i]<< (8-dec*2))+dec;
          i++;
       }
       pseq = realloc(seq,i+1);    // redimensionner le tableau
       if(pseq != NULL) seq = pseq;
       seq[i] = '\0';             // et le terminer
    }
    return seq;
}
//
 
int seqlen(unsigned char* seq){
    int lg = strlen(seq);
    int dec = seq[lg-1] & MASK_0;
 
    if (dec == 3) return lg*3;
    return (lg-1)*3+dec;
}
 
int seqcmp(unsigned char* seq1, int position1, unsigned char* seq2, int position2, int costMismatch){
 
    int offset2=position2%3, offset1 = position1%3;
    int nb1, nb2;
 
    switch (offset1){
        case 0: nb1 = seq1[position1/3] & MASK_3;
        case 1: nb1 = seq1[position1/3] & MASK_2;
        case 2: nb1 = seq1[position1/3] & MASK_1;
        default: return -1;
    }
 
    switch (offset2){
        case 0: nb2 = seq2[position2/3] & MASK_3;
        case 1: nb2 = seq2[position2/3] & MASK_2;
        case 2: nb2 = seq2[position2/3] & MASK_1;
        default: return -1;
    }
 
    if (offset2>offset1){
        if (nb2 == pow(3,offset2-offset1)*nb1) return 0;
        return costMismatch;
    }
    if (nb2 == pow(3,offset1-offset2)*nb1) return 0;
    return costMismatch;
}
 
unsigned char* reverseSequence(unsigned char* seq, int complement){
    int i, lg, dec;
    int seq_lg = seqlen(seq);
 
    //Calcul de la longueur du tableau de caractère
    if (seq_lg%3 == 0) lg = (seq_lg/3);
    else lg = (seq_lg/3)+1;
 
    dec = seq_lg%3;
    if (dec == 0) dec = 3;
 
    char *newSeq = calloc(lg+1, sizeof(unsigned char));
 
    if (newSeq != NULL){
        if (complement==0){
            for (i=0; i<lg; i++){
                newSeq[i] = ~(seq[i]) & MASK_1_2_3 | MASK_0;
            }
            newSeq[i-1] = (newSeq[i-1] & (MASK_0_1_2_3 << 8-dec*2)) | dec;
            newSeq[lg] = '\0';
            return newSeq;
        }
 
        for (i=0; i<lg/2; i++){
            newSeq[i] = reverseCompChar(seq[lg-i-1]);
            newSeq[lg-i-1] = reverseCompChar(seq[i]);
        }
 
        if (lg%2 == 1) newSeq[i] = reverseCompChar(seq[i]);
 
        newSeq[lg-1] = newSeq[lg-1] | MASK_0;
        newSeq[lg] = '\0';
 
        if (seq_lg%3 != 0){
            char* res = subSequence(newSeq, 3-(seq_lg%3), lg*3-1);
            free(newSeq);
            return res;
        }
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* subSequence(unsigned char* seq, int posDeb, int posFin){
 
    int i, decTmp = (posDeb%3)*2, k=0, lg, newSeq_lg, seq_lg, dec;
    unsigned char tmp1;
    unsigned char* newSeq;
 
    //Calcul de la longueur du tableau de caractère
    newSeq_lg = posFin-posDeb+1;
    seq_lg = seqlen(seq);
    dec = newSeq_lg%3;
 
    if (dec == 0) lg = newSeq_lg/3;
    else lg = newSeq_lg/3+1;
    newSeq = calloc(lg+1,sizeof(unsigned char));
 
    if (newSeq != NULL){
        for (i = posDeb/3; i<posFin/3; i++){
            tmp1 = (seq[i] & MASK_1_2_3) << decTmp;
            if (seq_lg > 3) tmp1 = tmp1 | (seq[i+1] >> (6-decTmp));
            newSeq[k] = tmp1 & MASK_1_2_3 | MASK_0;
            k++;
        }
 
        if (posDeb%3!=0){
            if (k>=lg){
                k--;
                newSeq[k] = (newSeq[k] & MASK_1_2_3) & (MASK_0_1_2_3 << 8-dec*2);
            }
            else newSeq[k] = ((seq[i] & MASK_1_2_3) << decTmp) & (MASK_0_1_2_3 << 8-dec*2);
        }
        else{
            if ((lg != 1) || (lg==1 && ((posDeb/3)==(posFin/3)))){
                if (dec == 0) newSeq[k] = seq[i] & MASK_1_2_3;
                else newSeq[k] = (seq[i] & MASK_1_2_3) & (MASK_0_1_2_3 << 8-dec*2);
            }
        }
 
        if (dec == 0) dec = 3;
        if (dec!=0) newSeq[lg-1] = newSeq[lg-1] | dec;
        newSeq[lg] = '\0';
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* concatSequence(unsigned char* seq1, unsigned char* seq2){
    int lgTab1, lgTab2, i, k, decTmp, lgTmp;
    int seq_lg, seq1_lg, seq2_lg;
    char* seq;
 
    // Oninitialise la nouvelle séquence
    seq1_lg = seqlen(seq1);
    seq2_lg = seqlen(seq2);
 
    seq_lg = seq1_lg + seq2_lg;
 
    if (seq_lg%3 == 0) lgTmp = seq_lg/3;
    else lgTmp = (seq_lg/3)+1;
 
    // Calcul de la longueur du tableau de la séquence 1
    if (seq1_lg%3 == 0) lgTab1 = seq1_lg/3;
    else lgTab1 = (seq1_lg/3)+1;
 
    // Calcul de la longueur du tableau de la séquence 2
    if (seq2_lg%3 == 0) lgTab2 = seq2_lg/3;
    else lgTab2 = (seq2_lg/3)+1;
 
    seq = calloc(lgTmp+1, sizeof(unsigned char));
 
    // Recopie du tableau de la séquence 1 dans la nouvelle séquence
    for (i=0; i<lgTab1; i++){
        seq[i] = seq1[i];
    }
 
    decTmp = (seq1[i-1] & MASK_0) * 2;
    seq[i-1] = seq[i-1] & MASK_1_2_3 /*| MASK_0*/;
 
    // Si aucun décalage dans la séquence 1, je recopie le tableau de la séquence 2 à la suite du 1 dans la nouvelle séquence
    if (decTmp == 6){
        seq[i-1] = seq[i-1] | MASK_0;
        for (k=0; k<lgTab2; k++){
            seq[i] = seq2[k];
            i++;
        }
    }
    else {
        // Gestion du décalage entre la dernière ligne du premier tableau et la première ligne du dernier tableau
 
        seq[i-1] = seq[i-1] | (seq2[0] >> decTmp) /*| MASK_0*/;
        if (i < lgTmp) seq[i] = (seq2[0] & MASK_1_2_3) << (6-decTmp);
 
        // Recopie du deuxième tableau avec gestion du décalage
        for (k=0; k<lgTab2-1; k++){
            seq[i] = seq[i] | (seq2[k+1] >> decTmp) /*| MASK_0*/;
            seq[i+1] = (seq2[k+1] & MASK_1_2_3) << (6-decTmp);
            i++;
        }
 
        if (seq_lg%3 == 0) seq[lgTmp-1] = seq[lgTmp-1] | MASK_0;
        else seq[lgTmp-1] = seq[lgTmp-1] | (seq_lg%3);
    }
 
    seq[lgTmp]= '\0';
 
    return seq;
}
 
unsigned char reverseCompChar(unsigned char c){
    int k;
    unsigned char tmp = 0;
 
    for (k=0; k<3; k++){
        switch(k){
            case 0: tmp = ((~(c >> 2) & 3) | tmp) << 2; break;
            case 1: tmp = ((~(c >> 4) & 3) | tmp) << 2; break;
            case 2: tmp = ((~(c >> 6) & 3) | tmp) << 2; break;
        }
    }
 
    tmp = tmp | MASK_0;
 
    return tmp;
}

Ne serait-il pas plus optimisé d'enregistrer 4 nucléotides par char ?

Tu n'aurais ainsi qu'à noter le nombre de nucléotides et/ou un pointeur vers le premier char non complet ainsi que le nombre de nucléotides modulo 4.

Ensuite tu peux allouer un tableau en 2D : conversion[nbNucléotideSurLeChar][Char];

Ce qui te donnerais :
Code:

1 2 3 4 for(char * p = adn.begin ; p != adn.fin ; ++p) printf("%s", conversion[4][*p]); if(adn.modulo) printf("%s", conversion[adn.modulo][*p]);
Pour l'opération inverse, je pense que le plus rapide est de jouer avec les pointeurs :
Code:

1 2 3 4 5 6 7 8 9 char possibilité [256]; char * dernierChoix[64][256] = {0}; char ** troisiemeChoix[16][256] = {0}; char *** deuxiemeChoix[4][256] = {0}; char **** premierChoix[256] = {0}; premierChoix['T'] = deuxièmeChoix[1]; //etc...
Ainsi au final tu pourras faire : char p = premierChoix['A']['T']['A']['G'];
Ce devrait être plus rapide que de faire des 'if'.

Certes pour l'utilisation de cette méthode, il faudra plus de mémoire mais pour le stockage de tes chaînes, tu devrais gagner 25% de place donc sur de grosses chaînes c'est négligeable.

03/08/2012, 14h24
Wizard50
Eh bien en théorie, 4 nucléotides sur un char c'est parfait, mais comme je l'ai dit précédemment, la séquence "AAAA" deviendra alors le char "\0", et comme dans de nombreux fonctions et structure (tel que les HashTables que j'utilise), on utilise des fonctions string tel que strlen, strcpy, etc ... tout buguera. Par exemple, un strlen() sur la séquence "AAAAGTCTGATGCTAAATGATC" renverra tout simplement 0.

De surcroit, il faut à coté du tableau indiquer soit sa longueur, soit le décalage de la dernière case. Indiquer le longueur du tableau a coté du tableau dans un structure, sa revient a faire
Code:

1 2 3 4 5 typedef struct { unsigned char* seq; int lg; } sequence;
Mais sur un OS 64 bits comme le mien, le processeur va ajouter 4 octets de bourrage pour aligner les données de la structure sur 8 bits (taille du bus de données). On peut aussi ajouter le décalage dans cette structure pour tout aligner mais au final, on perd de l'information quand même.

Sinon on peut utiliser un tableau 2D comme tu l'as dit mais hélas plusieurs désavantages :
- Sa ne passera pas comme clé d'une HashTable
- un octet pour le décalage + un octet pour stocker 4 nucléotides, sa reviens a 4 bits par nucléotides ce qui n'est pas vraiment efficace.

Merci quand même de tes propositions :)
03/08/2012, 17h20
Neckara
Citation:

Envoyé par Wizard50

mais comme je l'ai dit précédemment, la séquence "AAAA" deviendra alors le char "\0", et comme dans de nombreux fonctions et structure (tel que les HashTables que j'utilise), on utilise des fonctions string tel que strlen, strcpy, etc ... tout buguera. Par exemple, un strlen() sur la séquence "AAAAGTCTGATGCTAAATGATC" renverra tout simplement 0.

C'est pour cela que je t'ai donné un exemple de parcours.
il faudra juste redéfinir certaines fonctions de bases mais c'est très rapide.
Citation:
Envoyé par Wizard50

De surcroit, il faut à coté du tableau indiquer soit sa longueur, soit le décalage de la dernière case. Indiquer le longueur du tableau a coté du tableau dans un structure, sa revient a faire

Code:

1 2 3 4 5 typedef struct { unsigned char* seq; int lg; } sequence;

Mais sur un OS 64 bits comme le mien, le processeur va ajouter 4 octets de bourrage pour aligner les données de la structure sur 8 bits (taille du bus de données). On peut aussi ajouter le décalage dans cette structure pour tout aligner mais au final, on perd de l'information quand même.
Je ne comprend pas tout à fait ce que tu veux dire. Même si le processeur ajoute 4 octets de bourrage, c'est rien en comparaison des kio voir mio gagnés en mettant 4 nucléotides par char au lieu de 3.

Citation:

Envoyé par Wizard50

Sinon on peut utiliser un tableau 2D comme tu l'as dit mais hélas plusieurs désavantages :
- Sa ne passera pas comme clé d'une HashTable
- un octet pour le décalage + un octet pour stocker 4 nucléotides, sa reviens a 4 bits par nucléotides ce qui n'est pas vraiment efficace.

Je ne comprend pas tout. Le but du tableau 2D est de rendre beaucoup plus rapide les conversions entre données binaires et la chaîne ADN lisible par un humain. Sa taille est plus que négligeable en comparaison de la taille nécessaire pour stocker la chaîne d'ADN.

Mmmmmmmmh je pense avoir compris ton système. En soi la rapidité n'est pas un problème dans mon appli : que j'utilise un string standard ou un système pour coder 3 nucléotides sur un octet, le temps est généralement au pire (vraiment pire), à trois secondes près sur mes plus grosses données.

Mais si je comprend bien avec ton système, les HashTables que j'utilise doivent utiliser ma méthode de parcours et non pas une méthode de parcours de string standard. Le problème, c'est que l'on m'a demandé de concevoir ce système sans avoir à modifier la librarie de HashTable :(
Bien entendu, je comprend ta réflexion a propos des bits de padding et tu as raison.

De toute manière, un imprévu est arrivé : j'ai du prendre en compte le code d'ambiguité. C'est très simple : en plus des 4 nucléotides standard pour l'ADN et l'ARN (A, C, G et T/U), il existe 11 autres lettres qui designe une ambiguité, comme par exemple la lettre N qui désigne la nucléotide A ou C ou G ou T/U.
Cela me fait au total 15 symboles, ce qui est parfait : j'utilise 4 bits pour stocker une nucléotide, en utilisant les chiffres 1 à 15, et je stock deux nucléotides par octet. De cette manière, je ne peut pas avoir un octet avec '\0' et si les 4 bits de poids faible de mon dernier octet sont à 0, c'est que j'ai un décalage de 1 à la fin.

Le code s'en retrouve considérablement raccourci et plus simple.

Il y a juste que voilà, ma curiosité me titille et je n'ai pu m'enpecher de comparer mes deux systèmes (qui marche toujours très bien, pas de fuite mémoire) : le système qui code 3 nucléotides sur un octet utilise toujours au moins autant de mémoire que le système qui code 2 nucléotides sur un octet, et souvent plus même ! Hormis ces deux systèmes, le restant du code est le même, je n'y touche pas. J'ai essayé avec plusieurs taille de chaine et le résultat est toujours le même.
Pour mesurer la mémoire, je me suis basé sur le RSS de mon programme, qui d'après ce que j'ai compris, donne la taille en page virtuelle du code + la taille des données stockés dans la RAM.
Juste le fait d'utiliser les méthodes createSequence() et subSequence() suffit a me montrer que utiliser 2 nucléotides sur un octet est plus optimal en mémoire que 3 nucléotides sur un octet ...

Y'a t'il quelque chose qui m'échappe ? Fait-je quelque chose de travers ? Le RSS n'est pas la bonne mesure ?
Si quelqu'un pouvait me filer un tuyau, il me soulagerait de bien des tourments et je l'en remercie :)

Code pour stocker 2 nucléotides sur un octet :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "sequence_4bits.h"
#include <math.h>
 
#define MASK_0_1 15
#define MASK_1_0 240
#define MASK_1_1 255
 
enum {NUCLEO_0, NUCLEO_A, NUCLEO_C, NUCLEO_M, NUCLEO_R, NUCLEO_V, NUCLEO_H, NUCLEO_W, NUCLEO_S, NUCLEO_D, NUCLEO_B, NUCLEO_Y, NUCLEO_K, NUCLEO_G, NUCLEO_T_U, NUCLEO_N, NUCLEO_IGN};
//--------------------------------------------------------------
unsigned char* createSequence(char *str)
{
    int lg, i, dec;
    int nucleo;
    unsigned char * seq;
    unsigned char * pseq;
    //Initialisation du tableau de caractère
    lg = strlen(str);
    dec = lg%2;
    if (dec == 0) lg = lg/2;
    else lg = lg/2+1;
    seq = calloc(lg+1, sizeof(unsigned char));
    if (seq != NULL)
    {
      //lire la séquence
       for(i = 0, dec = 0 ;*str != '\0'; str++)
       {
          switch(*str)
          {
             case 'a' :
             case 'A' : nucleo = NUCLEO_A;   break;
             case 'c' :
             case 'C' : nucleo = NUCLEO_C;   break;
             case 'g' :
             case 'G' : nucleo = NUCLEO_G;   break;
             case 'u' :
             case 'U' :
             case 't' :
             case 'T' : nucleo = NUCLEO_T_U;   break;
             case 'm' :
             case 'M' : nucleo = NUCLEO_M;   break;
             case 'r' :
             case 'R' : nucleo = NUCLEO_R;   break;
             case 'w' :
             case 'W' : nucleo = NUCLEO_W;   break;
             case 's' :
             case 'S' : nucleo = NUCLEO_S;   break;
             case 'y' :
             case 'Y' : nucleo = NUCLEO_Y;   break;
             case 'k' :
             case 'K' : nucleo = NUCLEO_K;   break;
             case 'v' :
             case 'V' : nucleo = NUCLEO_V;   break;
             case 'h' :
             case 'H' : nucleo = NUCLEO_H;   break;
             case 'd' :
             case 'D' : nucleo = NUCLEO_D;   break;
             case 'b' :
             case 'B' : nucleo = NUCLEO_B;   break;
             case 'n' :
             case 'N' : nucleo = NUCLEO_N;   break;
             case ' ' :
             case '\n': nucleo = NUCLEO_IGN; break;
             default  : printf("Invalid character encountered : - Operation aborted");
                        free(seq);
                        exit(0);
          }
          if(nucleo != NUCLEO_IGN) // ne traiter que si c'est un caractère valide
          {
             if (dec==0){
                seq[i] = nucleo << 4;
                dec++;
             }
             else {
                seq[i] |=  nucleo;
                dec=0;
                i++;
             }
          }
       } // fin de lecture
       if (dec==1) i++;
       pseq = realloc(seq,i+1);    // redimensionner le tableau
       if(pseq != NULL) seq = pseq;
       seq[i] = '\0';             // et le terminer
    }
    return seq;
}
//
 
int seqlen(unsigned char* seq){
    int lg = strlen(seq);
    int dec = seq[lg-1] & MASK_0_1;
 
    if (dec == 0) return (lg-1)*2+1;
    return 2*lg;
}
 
int seqcmp(unsigned char* seq1, int position1, unsigned char* seq2, int position2, int costMismatch){
    int nb1, nb2;
 
    switch (position1%2){
        case 0: nb1 = seq1[position1/2] & MASK_0_1;
        case 1: nb1 = (seq1[position1/2] & MASK_1_0) >> 4;
        default: return -1;
    }
 
    switch (position2%2){
        case 0: nb2 = seq2[position2/2] & MASK_0_1;
        case 1: nb2 = (seq2[position2/2] & MASK_1_0) >> 4;
        default: return -1;
    }
 
    if (nb1 == nb2) return 1;
    return costMismatch;
}
 
unsigned char* reverseSequence(unsigned char* seq, int complement){
    int i, lg, seq_lg = seqlen(seq);
 
    //Calcul de la longueur du tableau de caractère
    lg = strlen(seq);
 
    char *newSeq = calloc(lg+1, sizeof(unsigned char));
 
    if (newSeq != NULL){
        if (complement==0){
            for (i=0; i<lg; i++){
                newSeq[i] = reverseCompChar(seq[i]);
            }
            if (seq_lg%2==1) newSeq[i-1] = newSeq[i-1] & MASK_1_0;
            newSeq[lg] = '\0';
            return newSeq;
        }
 
        if (lg > 1){
            for (i=0; i<lg/2; i++){
                newSeq[i] = reverseCompChar(seq[lg-i-1]);
                newSeq[lg-i-1] = reverseCompChar(seq[i]);
            }
 
            if (lg%2 == 1) newSeq[i] = reverseCompChar(seq[i]);
            newSeq[lg] = '\0';
 
            if (seq_lg%2 == 1){
                char* res = subSequence(newSeq, 1, (lg*2)-1);
                free(newSeq);
                return res;
            }
        }
        else{
            newSeq[i] = reverseCompChar(seq[i]);
            if (seq_lg%2==1) newSeq[i] = newSeq[i] << 4;
            newSeq[lg] = '\0';
        }
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* subSequence(unsigned char* seq, int posDeb, int posFin){
 
    int i, k=0, lg, newSeq_lg, seq_lg, dec;
    unsigned char tmp1;
    unsigned char* newSeq;
 
    //Calcul de la longueur du tableau de caractère
    newSeq_lg = posFin-posDeb+1;
    seq_lg = seqlen(seq);
    dec = newSeq_lg%2;
 
    if (dec == 0) lg = newSeq_lg/2;
    else lg = newSeq_lg/2+1;
    newSeq = calloc(lg+1,sizeof(unsigned char));
 
    if (newSeq != NULL){
        if (posDeb%2 == 0){
            for (i = posDeb/2; i<posFin/2; i++){
                newSeq[k] = seq[i];
                k++;
            }
            if (dec == 1) newSeq[k] = seq[i] & MASK_1_0;
            else newSeq[k] = seq[i];
        }
        else {
            for (i = posDeb/2; i<posFin/2; i++){
                tmp1 = seq[i] << 4;
                if (seq_lg > 2) tmp1 = tmp1 | (seq[i+1] >> 4);
                newSeq[k] = tmp1;
                k++;
            }
            if (dec == 1) newSeq[k] = seq[i] << 4;
            //else newSeq[k] = seq[i];
        }
 
        newSeq[lg] = '\0';
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* concatSequence(unsigned char* seq1, unsigned char* seq2){
    int lgTab1, lgTab2, i, k, lgTmp, decTmp = 0;
    int seq_lg;
    unsigned char* seq;
 
    seq_lg = seqlen(seq1) + seqlen(seq2);
 
    if (seq_lg%2 == 0) lgTmp = seq_lg/2;
    else lgTmp = (seq_lg/2)+1;
 
    // Calcul de la longueur du tableau de la séquence 1
    lgTab1 = strlen(seq1);
    lgTab2 = strlen(seq2);
 
    seq = calloc(lgTmp+1, sizeof(unsigned char));
 
    // Recopie du tableau de la séquence 1 dans la nouvelle séquence
    for (i=0; i<lgTab1; i++){
        seq[i] = seq1[i];
    }
 
    if ((seq1[i-1] & MASK_0_1) == 0) decTmp = 1;
 
    // Si aucun décalage dans la séquence 1, je recopie le tableau de la séquence 2 à la suite du 1 dans la nouvelle séquence
    if (decTmp == 0){
        for (k=0; k<lgTab2; k++){
            seq[i] = seq2[k];
            i++;
        }
    }
    else {
        // Gestion du décalage entre la dernière ligne du premier tableau et la première ligne du dernier tableau
 
        seq[i-1] = seq[i-1] | (seq2[0] >> 4);
        if (i < lgTmp) seq[i] = seq2[0] << 4;
 
        // Recopie du deuxième tableau avec gestion du décalage
        for (k=0; k<lgTab2-1; k++){
            seq[i] = seq[i] | (seq2[k+1] >> 4);
            seq[i+1] = seq2[k+1] << 4;
            i++;
        }
    }
 
    seq[lgTmp]= '\0';
 
    return seq;
}
 
unsigned char reverseCompChar(unsigned char c){
    int k;
    unsigned char tmp = 0;
 
    for (k=0; k<2; k++){
        switch(k){
            case 0: tmp = reverseNucl(c) << 4; break;
            case 1: tmp = reverseNucl(c >> 4) | tmp; break;
        }
    }
 
    return tmp;
}
 
unsigned char reverseNucl(unsigned char c){
    switch (c){
        case 7:
        case 8:
        case 15: return c;
        default: return ~c & MASK_0_1;
    }
}

Code pour stocker 3 nucléotides sur un octet :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
 
 
#include "sequence.h"
#include <math.h>
 
#define MASK_0 3
#define MASK_1 12
#define MASK_2 48
#define MASK_3 192
#define MASK_3_2 240
#define MASK_1_2_3 252
#define MASK_0_1_2_3 255
 
enum {NUCLEO_A, NUCLEO_C, NUCLEO_G, NUCLEO_T, NUCLEO_IGN};
//--------------------------------------------------------------
unsigned char* createSequence(char *str)
{
    int lg, lg_str, i=0, j=0, dec;
    int nucleo;
    unsigned char * seq;
    unsigned char * pseq;
    //Initialisation du tableau de caractère
    lg_str = strlen(str);
    dec = lg_str%3;
    if (dec == 0) lg = lg_str/3;
    else {
        lg = lg_str/3+1;
        dec = 0;
    }
 
    seq = calloc(lg+1, sizeof(unsigned char));
    if (seq != NULL)
    {
      //lire la séquence
       while (j < lg_str)
       {
          switch(str[j])
          {
             case 'a' :
             case 'A' : nucleo = NUCLEO_A;   break;
             case 'c' :
             case 'C' : nucleo = NUCLEO_C;   break;
             case 'g' :
             case 'G' : nucleo = NUCLEO_G;   break;
             case 't' :
             case 'T' : nucleo = NUCLEO_T;   break;
             case ' ' :
             case '\n': nucleo = NUCLEO_IGN; break;
             default  : printf("Invalid character encountered : - Operation aborted");
                        free(seq);
                        exit(0);
          }
          if(nucleo != NUCLEO_IGN) // ne traiter que si ACGT
          {
             seq[i] = (seq[i]<<2) | nucleo;
             dec++;
             if(dec==3) // octet complet
             {
                seq[i] = (seq[i]<<2)+3;
                i++;
                dec= 0;
             }
          }
          j++;
       } // fin de lecture
       if(dec != 0)
       {
          seq[i] = (seq[i]<< (8-dec*2))+dec;
          i++;
       }
       //printf("i = %d, lg = %d \n", i, lg);
       if (i>lg){
           pseq = (unsigned char*)realloc(seq,(i+1)*sizeof(unsigned char));
           if(pseq != NULL) seq = pseq;
       }
       seq[i] = '\0';
    }
    return seq;
}
 
int seqlen(unsigned char* seq){
    int lg = strlen(seq);
    int dec = seq[lg-1] & MASK_0;
 
    if (dec == 3) return lg*3;
    return (lg-1)*3+dec;
}
 
int seqcmp(unsigned char* seq1, int position1, unsigned char* seq2, int position2, int costMismatch){
 
    int offset2=position2%3, offset1 = position1%3;
    int nb1, nb2;
 
    switch (offset1){
        case 0: nb1 = seq1[position1/3] & MASK_3;
        case 1: nb1 = seq1[position1/3] & MASK_2;
        case 2: nb1 = seq1[position1/3] & MASK_1;
        default: return -1;
    }
 
    switch (offset2){
        case 0: nb2 = seq2[position2/3] & MASK_3;
        case 1: nb2 = seq2[position2/3] & MASK_2;
        case 2: nb2 = seq2[position2/3] & MASK_1;
        default: return -1;
    }
 
    if (offset2>offset1){
        if (nb2 == pow(3,offset2-offset1)*nb1) return 0;
        return costMismatch;
    }
    if (nb2 == pow(3,offset1-offset2)*nb1) return 0;
    return costMismatch;
}
 
unsigned char* reverseSequence(unsigned char* seq, int complement){
    int i, lg, dec;
    int seq_lg = seqlen(seq);
 
    //Calcul de la longueur du tableau de caractère
    lg = strlen(seq);
    dec = seq_lg%3;
 
    char *newSeq = calloc(lg+1, sizeof(unsigned char));
 
    if (newSeq != NULL){
        if (complement==0){
            if (dec == 0) dec = 3;
            for (i=0; i<lg; i++){
                newSeq[i] = ~(seq[i]) & MASK_1_2_3 | MASK_0;
            }
            newSeq[i-1] = (newSeq[i-1] & (MASK_0_1_2_3 << 8-dec*2)) | dec;
            newSeq[lg] = '\0';
            return newSeq;
        }
 
        if (lg == 1){
            newSeq[0] = reverseCompChar(seq[lg-i-1]) | MASK_0;
            newSeq[1] = '\0';
        }
        else{
            for (i=0; i<lg/2; i++){
                newSeq[i] = reverseCompChar(seq[lg-i-1]);
                newSeq[lg-i-1] = reverseCompChar(seq[i]);
            }
 
            if (lg%2 == 1) newSeq[i] = reverseCompChar(seq[i]);
 
            newSeq[lg-1] = newSeq[lg-1] | MASK_0;
            newSeq[lg] = '\0';
        }
 
        if (dec != 0){
            unsigned char* res = subSequence(newSeq, 3-dec, (lg*3)-1);
            free(newSeq);
            return res;
        }
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* subSequence(unsigned char* seq, int posDeb, int posFin){
 
    int i, k=0, lg, newSeq_lg, dec;
    unsigned char* newSeq;
 
    //Calcul de la longueur du tableau de caractère
    newSeq_lg = posFin-posDeb+1;
    dec = newSeq_lg%3;
 
    if (dec == 0) lg = newSeq_lg/3;
    else lg = newSeq_lg/3+1;
    newSeq = calloc(lg+1,sizeof(unsigned char));
 
    if (newSeq != NULL){
        if (posDeb%3 == 0){
            for (i = posDeb/3; i<posFin/3; i++){
                newSeq[k] = seq[i];
                k++;
            }
            if (dec == 1) newSeq[k] = (seq[i] & MASK_3) + 1;
            else if (dec == 2) newSeq[k] = (seq[i] & MASK_3_2) + 2;
            else newSeq[k] = (seq[i] & MASK_1_2_3) + 3;
        }
        else {
            int decTmp = (posDeb%3)*2;
            int tmp1;
            for (i = posDeb/3; i<posFin/3; i++){
                tmp1 = (seq[i] & MASK_1_2_3) << decTmp;
                if (lg > 1) tmp1 = tmp1 | (seq[i+1] >> (6-decTmp));
                newSeq[k] = (tmp1 & MASK_1_2_3) | MASK_0;
                k++;
            }
            if (dec == 1) newSeq[k] = ((seq[i] & MASK_1_2_3) << 2) + 1;
            else if (dec == 2) newSeq[k] = ((seq[i] & MASK_1_2_3) << 4) + 2;
            else newSeq[k] = (seq[i] & MASK_1_2_3) + 3;
        }
 
        newSeq[lg] = '\0';
 
        return newSeq;
    }
    else {
        fprintf (stderr, "Insufficient memory\n");
        exit (EXIT_FAILURE);
    }
}
 
unsigned char* concatSequence(unsigned char* seq1, unsigned char* seq2){
    int lgTab1, lgTab2, i, k, decTmp, lgTmp;
    int seq_lg;
    unsigned char* seq;
 
    seq_lg = seqlen(seq1) + seqlen(seq2);
 
    if (seq_lg%3 == 0) lgTmp = seq_lg/3;
    else lgTmp = (seq_lg/3)+1;
 
    // Calcul de la longueur du tableau de la séquence 1
    lgTab1 = strlen(seq1);
    lgTab2 = strlen(seq2);
 
    seq = calloc(lgTmp+1, sizeof(unsigned char));
 
    // Recopie du tableau de la séquence 1 dans la nouvelle séquence
    for (i=0; i<lgTab1; i++){
        seq[i] = seq1[i];
    }
 
    decTmp = (seq1[i-1] & MASK_0) * 2;
    seq[i-1] = seq[i-1] & MASK_1_2_3 /*| MASK_0*/;
 
    // Si aucun décalage dans la séquence 1, je recopie le tableau de la séquence 2 à la suite du 1 dans la nouvelle séquence
    if (decTmp == 6){
        seq[i-1] = seq[i-1] | MASK_0;
        for (k=0; k<lgTab2; k++){
            seq[i] = seq2[k];
            i++;
        }
    }
    else {
        // Gestion du décalage entre la dernière ligne du premier tableau et la première ligne du dernier tableau
 
        seq[i-1] = seq[i-1] | (seq2[0] >> decTmp) /*| MASK_0*/;
        if (i < lgTmp) seq[i] = (seq2[0] & MASK_1_2_3) << (6-decTmp);
 
        // Recopie du deuxième tableau avec gestion du décalage
        for (k=0; k<lgTab2-1; k++){
            seq[i] = seq[i] | (seq2[k+1] >> decTmp) /*| MASK_0*/;
            seq[i+1] = (seq2[k+1] & MASK_1_2_3) << (6-decTmp);
            i++;
        }
 
        if (seq_lg%3 == 0) seq[lgTmp-1] = seq[lgTmp-1] | MASK_0;
        else seq[lgTmp-1] = seq[lgTmp-1] | (seq_lg%3);
    }
 
    seq[lgTmp]= '\0';
 
    return seq;
}
 
unsigned char reverseCompChar(unsigned char c){
    int k;
    unsigned char tmp = 0;
 
    for (k=0; k<3; k++){
        switch(k){
            case 0: tmp = ((~(c >> 2) & 3) | tmp) << 2; break;
            case 1: tmp = ((~(c >> 4) & 3) | tmp) << 2; break;
            case 2: tmp = ((~(c >> 6) & 3) | tmp) << 2; break;
        }
    }
 
    tmp = tmp | MASK_0;
 
    return tmp;
}

07/08/2012, 17h52
Neckara
Si tu es sous Linux, tu peux aussi utiliser la commande time :
http://man.developpez.com/man1/time.1.php

( %M pour la mémoire)

Sinon je suis très étonné par les différences de performances entre les deux codes.
Je pense tout de même qu'on peut optimiser un peu :
Code:

1 2 if (dec == 0) lg = lg/2; else lg = lg/2+1;
Code:

lg = lg/2 + (dec != 0)

Une autre chose, calloc alloue un espace en mémoire et initialise chaque octet à 0, il est donc beaucoup plus lent qu'un malloc, as-tu bien besoin de l'utiliser?
Code:

1 2 seq[i] = seq[i] | (seq2[k+1] >> 4); seq[i] |= (seq2[k+1] >> 4);//plus rapide car une copie en moins
Code:

1 2 3 4 5 6 for (k=0; k<2; k++){ switch(k){ case 0: tmp = reverseNucl(c) << 4; break; case 1: tmp = reverseNucl(c >> 4) | tmp; break; } }
hein??????
07/08/2012, 21h20
Wizard50

Merci pour toutes tes remarques, je vais essayer de les appliquer sur mon programme dès ce soir et vous redonner les résultats de la commande time au passage.

Pour ce qui est de la différence de performance entre les deux systèmes, je t'assure que il n'y a presque rien mais dès le prochain test que je fais, je te donne les résultats.

En ce qui concerne le dernier quote, je ne sais pas quoi te dire ... ou plutot comment me justifier si ce n'est beaucoup de fatigue en ce moment :)
08/08/2012, 11h31
Wizard50

Donc voilà, petit test avec un fichier de 2 millions de caractères.

Système codant 3 nucléotides sur 1 octet :
- Commande time avec %M : 640 512 Ko
- Commande smem a propos du RSS : 156,4 Mo
- Temps d'execution : 3,17 s

Système codant 2 nucléotides sur 1 octet :
- Commande time avec %M : 547 696 Ko
- Commande smem a propos du RSS : 133,7 Mo
- Temps d'execution : 3,07 s

Je n'y comprend plus rien ...