source: tags/arb_5.0/ARBDB/adseqcompr.c

Last change on this file was 6143, checked in by westram, 16 years ago
  • backport [6141] (parts changing code, but only strings and comments)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.8 KB
Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4
5#include <adlocal.h>
6#include <arbdbt.h>
7
8/* -------------------------------------------------------------------------------- */
9
10#define MAX_SEQUENCE_PER_MASTER 50 /* was 18 till May 2008 */
11
12#if defined(DEBUG)
13/* don't do optimize, only create tree and save to DB */
14/* #define SAVE_COMPRESSION_TREE_TO_DB */
15#endif /* DEBUG */
16
17/* -------------------------------------------------------------------------------- */
18
19typedef struct gb_seq_compr_tree {
20#ifdef FAKE_VTAB_PTR
21    virtualTable *dummy_virtual; /* simulate pointer to virtual-table used in AP_tree */
22#endif
23    GBT_TREE_ELEMENTS(struct gb_seq_compr_tree);
24
25    int index;                  /* master(inner nodes) or sequence(leaf nodes) index */
26    int sons;                   /* sons with sequence or masters (in subtree) */
27} GB_CTREE;
28
29typedef struct {
30    int            len;
31    char           used[256];
32    unsigned char *con[256];
33} GB_Consensus;
34
35typedef struct {
36    GBDATA *gbd;
37    int     master;
38} GB_Sequence;
39
40typedef struct {
41    GBDATA *gbd;
42    int     master;
43} GB_Master;
44
45/* -------------------------------------------------------------------------------- */
46
47GB_Consensus *g_b_new_Consensus(long len){
48    GB_Consensus *gcon = (GB_Consensus *)GB_calloc(sizeof(GB_Consensus),1);
49    int i;
50    unsigned char *data = (unsigned char *)GB_calloc(sizeof(char)*256,len);
51    gcon->len = len;
52
53    for (i=0;i<256;i++){
54        gcon->con[i] = data + len*i;
55    }
56    return gcon;
57}
58
59
60void g_b_delete_Consensus(GB_Consensus *gcon){
61    free((char *)gcon->con[0]);
62    free((char *)gcon);
63}
64
65
66void g_b_Consensus_add(GB_Consensus *gcon, unsigned char *seq, long seq_len){
67    int            i;
68    int            li;
69    int            c;
70    unsigned char *s;
71    int            last;
72    unsigned char *p;
73    int            eq_count;
74    const int      max_priority = 255/MAX_SEQUENCE_PER_MASTER; /* No overflow possible */
75
76    ad_assert(max_priority >= 1);
77
78    if (seq_len > gcon->len) seq_len = gcon->len;
79
80    /* Search for runs  */
81    s = seq;
82    last = 0;
83    eq_count = 0;
84    for (li = i = 0; i < seq_len; i++ ){
85        c = *(s++);
86        if (c == last){
87            continue;
88        }else{
89        inc_hits:
90            eq_count = i-li;
91            gcon->used[c] = 1;
92            p = gcon->con[last];
93            last = c;
94            if (eq_count <= GB_RUNLENGTH_SIZE) {
95                c = max_priority;
96                while (li < i) p[li++] += c;
97            }else{
98                c = max_priority * ( GB_RUNLENGTH_SIZE ) / eq_count;
99                if (c){
100                    while (li < i) p[li++] += c;
101                }else{
102                    while (li < i) p[li++] |= 1;
103                }
104            }
105        }
106    }
107    if (li<seq_len){
108        c = last;
109        i = seq_len;
110        goto inc_hits;
111    }
112}
113
114char *g_b_Consensus_get_sequence(GB_Consensus *gcon){
115    int pos;
116    unsigned char *s;
117    unsigned char *max = (unsigned char *)GB_calloc(sizeof(char), gcon->len);
118    int c;
119    char *seq = (char *)GB_calloc(sizeof(char),gcon->len+1);
120
121    memset(seq,'@',gcon->len);
122
123    for (c = 1; c<256;c++){ /* Find maximum frequency of non run */
124        if (!gcon->used[c]) continue;
125        s = gcon->con[c];
126        for (pos= 0;pos<gcon->len;pos++){
127            if (*s > max[pos]) {
128                max[pos] = *s;
129                seq[pos] = c;
130            }
131            s++;
132        }
133    }
134    free((char *)max);
135    return seq;
136}
137
138
139int g_b_count_leafs(GB_CTREE *node){
140    if (node->is_leaf) return 1;
141    node->gb_node = 0;
142    return (g_b_count_leafs(node->leftson) + g_b_count_leafs(node->rightson));
143}
144
145void g_b_put_sequences_in_container(GB_CTREE *ctree,GB_Sequence *seqs,GB_Master **masters,GB_Consensus *gcon){
146    if (ctree->is_leaf){
147        if (ctree->index >= 0) {
148            GB_CSTR data = GB_read_char_pntr(seqs[ctree->index].gbd);
149            long    len  = GB_read_string_count(seqs[ctree->index].gbd);
150            g_b_Consensus_add(gcon,(unsigned char *)data,len);
151        }
152    }
153    else if (ctree->index<0){
154        g_b_put_sequences_in_container(ctree->leftson,seqs,masters,gcon);
155        g_b_put_sequences_in_container(ctree->rightson,seqs,masters,gcon);
156    }
157    else {
158        GB_CSTR data = GB_read_char_pntr(masters[ctree->index]->gbd);
159        long    len  = GB_read_string_count(masters[ctree->index]->gbd);
160        g_b_Consensus_add(gcon,(unsigned char *)data,len);
161    }
162}
163
164void g_b_create_master(GB_CTREE *node, GB_Sequence *seqs, GB_Master **masters,
165                       int masterCount, int *builtMasters, int my_master,
166                       const char *ali_name, long seq_len, int *aborted)
167{
168    if (*aborted) {
169        return;
170    }
171
172    if (node->is_leaf){
173        if (node->index >= 0) {
174            GBDATA *gb_data = GBT_read_sequence(node->gb_node, ali_name);
175
176            seqs[node->index].gbd    = gb_data;
177            seqs[node->index].master = my_master;
178        }
179    }
180    else {
181        if (node->index>=0){
182            masters[node->index]->master = my_master;
183            my_master = node->index;
184        }
185        g_b_create_master(node->leftson,seqs,masters,masterCount,builtMasters,my_master,ali_name,seq_len, aborted);
186        g_b_create_master(node->rightson,seqs,masters,masterCount,builtMasters,my_master,ali_name,seq_len, aborted);
187        if (node->index>=0 && !*aborted) { /* build me */
188            char         *data;
189            GB_Consensus *gcon = g_b_new_Consensus(seq_len);
190
191            g_b_put_sequences_in_container(node->leftson,seqs,masters,gcon);
192            g_b_put_sequences_in_container(node->rightson,seqs,masters,gcon);
193
194            data = g_b_Consensus_get_sequence(gcon);
195
196            GB_write_string(masters[node->index]->gbd,data);
197            GB_write_security_write(masters[node->index]->gbd,7);
198
199            g_b_delete_Consensus(gcon);
200            free(data);
201
202            (*builtMasters)++;
203            *aborted |= GB_status(*builtMasters/(double)masterCount);
204        }
205    }
206}
207
208/* ------------------------------------- */
209/*      distribute master sequences      */
210/* ------------------------------------- */
211
212static void subtract_sons_from_tree(GB_CTREE *node, int subtract) {
213    while (node) {
214        node->sons -= subtract;
215        node        = node->father;
216    }
217}
218
219static int set_masters_with_sons(GB_CTREE *node, int wantedSons, int *mcount) {
220    if (!node->is_leaf) {
221        if (node->sons == wantedSons) {
222            /* insert new master */
223            ad_assert(node->index == -1);
224            node->index = *mcount;
225            (*mcount)++;
226
227            subtract_sons_from_tree(node->father, node->sons-1);
228            node->sons = 1;
229        }
230        else if (node->sons>wantedSons) {
231            int lMax = set_masters_with_sons(node->leftson,  wantedSons, mcount);
232            int rMax = set_masters_with_sons(node->rightson, wantedSons, mcount);
233
234            int maxSons = lMax<rMax ? rMax : lMax;
235            if (node->sons <= MAX_SEQUENCE_PER_MASTER && node->sons>maxSons) {
236                maxSons = node->sons;
237            }
238            return maxSons;
239        }
240    }
241    return node->sons <= MAX_SEQUENCE_PER_MASTER ? node->sons : 0;
242}
243
244static int maxCompressionSteps(GB_CTREE *node) {
245    if (node->is_leaf) {
246        return 0;
247    }
248
249    int left  = maxCompressionSteps(node->leftson);
250    int right = maxCompressionSteps(node->rightson);
251
252#if defined(SAVE_COMPRESSION_TREE_TO_DB)
253    freeset(node->name, 0);
254    if (node->index2 != -1) {
255        node->name = GBS_global_string_copy("master_%03i", node->index2);
256    }
257#endif /* SAVE_COMPRESSION_TREE_TO_DB */
258
259    return (left>right ? left : right) + (node->index == -1 ? 0 : 1);
260}
261
262static int init_indices_and_count_sons(GB_CTREE *node, int *scount, const char *ali_name) {
263    if (node->is_leaf){
264        if (node->gb_node ==0 || !GBT_read_sequence(node->gb_node,(char *)ali_name)) {
265            node->index = -1;
266            node->sons  = 0;
267        }
268        else {
269            node->index = *scount;
270            node->sons  = 1;
271            (*scount)++;
272        }
273    }
274    else {
275        node->index = -1;
276        node->sons  =
277            init_indices_and_count_sons(node->leftson, scount, ali_name) +
278            init_indices_and_count_sons(node->rightson, scount, ali_name);
279    }
280    return node->sons;
281}
282
283static void distribute_masters(GB_CTREE *tree, int *mcount, int *max_masters) {
284    int wantedSons = MAX_SEQUENCE_PER_MASTER;
285    while (wantedSons >= 2) {
286        int maxSons = set_masters_with_sons(tree, wantedSons, mcount);
287        wantedSons = maxSons;
288    }
289    ad_assert(tree->sons == 1);
290
291    ad_assert(tree->index != -1);
292    *max_masters = maxCompressionSteps(tree);
293}
294
295/* -------------------------------------------------------------------------------- */
296
297static GB_INLINE long g_b_read_number2(const unsigned char **s) {
298    unsigned int c0,c1,c2,c3,c4;
299    c0 = (*((*s)++));
300    if (c0 & 0x80){
301        c1 = (*((*s)++));
302        if (c0 & 0x40) {
303            c2 = (*((*s)++));
304            if (c0 & 0x20) {
305                c3 = (*((*s)++));
306                if (c0 &0x10) {
307                    c4 = (*((*s)++));
308                    return c4 | (c3<<8) | (c2<<16) | (c1<<8);
309                }else{
310                    return (c3) | (c2<<8 ) | (c1<<16) | ((c0 & 0x0f)<<24);
311                }
312            }else{
313                return (c2) | (c1<<8) | ((c0 & 0x1f)<<16);
314            }
315        }else{
316            return (c1) | ((c0 & 0x3f)<<8);
317        }
318    }else{
319        return c0;
320    }
321}
322
323static GB_INLINE void g_b_put_number2(int i, unsigned char **s) {
324    int j;
325    if (i< 0x80 ){ *((*s)++)=i;return; }
326    if (i<0x4000) {
327        j = (i>>8) | 0x80;
328        *((*s)++)=j;
329        *((*s)++)=i;
330    }else if (i<0x200000) {
331        j = (i>>16) | 0xC0;
332        *((*s)++)=j;
333        j = (i>>8);
334        *((*s)++)=j;
335        *((*s)++)=i;
336    } else if (i<0x10000000) {
337        j = (i>>24) | 0xE0;
338        *((*s)++)=j;
339        j = (i>>16);
340        *((*s)++)=j;
341        j = (i>>8);
342        *((*s)++)=j;
343        *((*s)++)=i;
344    }
345}
346
347
348char *gb_compress_seq_by_master(const char *master,int master_len,int master_index,
349                                GBQUARK q, const char *seq, long seq_len,
350                                long *memsize, int old_flag){
351    unsigned char *buffer;
352    int            rest = 0;
353    unsigned char *d;
354    int            i,cs,cm;
355    int            last;
356    long           len  = seq_len;
357
358    d = buffer = (unsigned char *)GB_give_other_buffer(seq,seq_len);
359
360    if (seq_len > master_len){
361        rest = seq_len - master_len;
362        len = master_len;
363    }
364
365    last = -1000;       /* Convert Sequence relative to Master */
366    for( i = len; i>0; i--){
367        cm = *(master++);
368        cs = *(seq++);
369        if (cm==cs && cs != last){
370            *(d++) = 0;
371            last = 1000;
372        }else{
373            *(d++) = cs;
374            last = cs;
375        }
376    }
377    for( i = rest; i>0; i--){
378        *(d++) = *(seq++);
379    }
380
381    {               /* Append run length compression method */
382        unsigned char *buffer2;
383        unsigned char *dest2;
384        buffer2 = dest2 = (unsigned char *)GB_give_other_buffer((char *)buffer,seq_len+100);
385        *(dest2++) = GB_COMPRESSION_SEQUENCE | old_flag;
386
387        g_b_put_number2(master_index,&dest2); /* Tags */
388        g_b_put_number2(q,&dest2);
389
390        gb_compress_equal_bytes_2((char *)buffer,seq_len,memsize,(char *)dest2); /* append runlength compressed sequences to tags */
391
392        *memsize = *memsize + (dest2-buffer2);
393        return (char *)buffer2;
394    }
395}
396
397char *gb_compress_sequence_by_master(GBDATA *gbd,const char *master,int master_len,int master_index,
398                                     GBQUARK q, const char *seq, int seq_len, long *memsize)
399{
400    long  size;
401    char *is  = gb_compress_seq_by_master(master,master_len,master_index,q,seq,seq_len,&size,GB_COMPRESSION_LAST);
402    char *res = gb_compress_data(gbd,0,is,size,memsize, ~(GB_COMPRESSION_DICTIONARY|GB_COMPRESSION_SORTBYTES|GB_COMPRESSION_RUNLENGTH),GB_TRUE);
403    return res;
404}
405
406static GB_ERROR compress_sequence_tree(GBDATA *gb_main, GB_CTREE *tree, const char *ali_name){
407    GB_ERROR error      = 0;
408    long     ali_len    = GBT_get_alignment_len(gb_main, ali_name);
409    int      main_clock = GB_read_clock(gb_main);
410
411    if (ali_len<0){
412        GB_ERROR warning = GBS_global_string("Skipping alignment '%s' (not a valid alignment; len=%li).", ali_name, ali_len);
413        GB_information(warning);
414    }
415    else {
416        int leafcount = g_b_count_leafs(tree);
417        if (!leafcount) {
418            error = GB_export_errorf("Tree contains no sequences with data in '%s'", ali_name);
419        }
420        else {
421            /* Distribute masters in tree */
422            int mastercount   = 0;
423            int max_compSteps = 0; // in one branch
424            int seqcount      = 0;
425
426            GB_status2("Create master sequences");
427            GB_status(0);
428
429            init_indices_and_count_sons(tree, &seqcount, ali_name);
430            distribute_masters(tree, &mastercount, &max_compSteps);
431
432#if defined(SAVE_COMPRESSION_TREE_TO_DB)
433            {
434                error = GBT_link_tree((GBT_TREE*)tree, gb_main, 0, NULL, NULL);
435                if (!error) error = GBT_write_tree(gb_main,0,"tree_compression_new",(GBT_TREE*)tree);
436                GB_information("Only generated compression tree (do NOT save DB anymore)");
437                return error;
438            }
439#endif /* SAVE_COMPRESSION_TREE_TO_DB */
440
441            // detect degenerated trees
442            {
443                int min_masters   = ((seqcount-1)/MAX_SEQUENCE_PER_MASTER)+1;
444                int min_compSteps = 1;
445                {
446                    int m = min_masters;
447                    while (m>1) {
448                        m            = ((m-1)/MAX_SEQUENCE_PER_MASTER)+1;
449                        min_masters += m;
450                        min_compSteps++;
451                    }
452                }
453
454                int acceptable_masters   = (3*min_masters)/2; // accept 50% overhead
455                int acceptable_compSteps = 11*min_compSteps; // accept 1000% overhead
456
457                if (mastercount>acceptable_masters || max_compSteps>acceptable_compSteps) {
458                    GB_warningf("Tree is ill-suited for compression (cause of deep branches)\n"
459                                "                    Used tree   Optimal tree   Overhead\n"
460                                "Compression steps       %5i          %5i      %4i%% (speed)\n"
461                                "Master sequences        %5i          %5i      %4i%% (size)\n"
462                                "If you like to restart with a better tree,\n"
463                                "press 'Abort' to stop compression",
464                                max_compSteps, min_compSteps, (100*max_compSteps)/min_compSteps-100,
465                                mastercount, min_masters, (100*mastercount)/min_masters-100);
466                }
467            }
468
469            ad_assert(mastercount>0);
470
471            {
472                GBDATA              *gb_master_ali     = 0;
473                GBDATA              *old_gb_master_ali = 0;
474                GB_Sequence         *seqs              = 0;
475                GB_MAIN_TYPE        *Main              = GB_MAIN(gb_main);
476                GBQUARK              ali_quark         = gb_key_2_quark(Main,ali_name);
477                unsigned long long   sumorg            = 0;
478                unsigned long long   sumold            = 0;
479                unsigned long long   sumnew            = 0;
480                GB_Master          **masters           = (GB_Master **)GB_calloc(sizeof(*masters),leafcount);
481                int                  si;
482
483                {
484                    char *masterfoldername = GBS_global_string_copy("%s/@master_data/@%s",GB_SYSTEM_FOLDER,ali_name);
485                    old_gb_master_ali      = GB_search(gb_main, masterfoldername,GB_FIND);
486                    free(masterfoldername);
487                }
488
489                // create masters
490                if (!error) {
491                    {
492                        char   *master_data_name = GBS_global_string_copy("%s/@master_data",GB_SYSTEM_FOLDER);
493                        char   *master_name      = GBS_global_string_copy("@%s",ali_name);
494                        GBDATA *gb_master_data   = gb_search(gb_main, master_data_name,GB_CREATE_CONTAINER,1);
495
496                        /* create a master container, the old is deleted as soon as all sequences are compressed by the new method*/
497                        gb_master_ali = gb_create_container(gb_master_data, master_name);
498                        GB_write_security_delete(gb_master_ali,7);
499
500                        free(master_name);
501                        free(master_data_name);
502                    }
503                    for (si = 0; si<mastercount; si++) {
504                        masters[si]      = (GB_Master *)GB_calloc(sizeof(GB_Master),1);
505                        masters[si]->gbd = gb_create(gb_master_ali,"@master",GB_STRING);
506                    }
507                    seqs = (GB_Sequence *)GB_calloc(sizeof(*seqs),leafcount);
508                    {
509                        int builtMasters = 0;
510                        int aborted      = 0;
511                        GB_status2("Building %i master sequences", mastercount);
512                        g_b_create_master(tree,seqs,masters,mastercount,&builtMasters,-1,ali_name,ali_len, &aborted);
513                        if (aborted) error = "User abort";
514                    }
515                }
516
517                /* Compress sequences in tree */
518                if (!error) {
519                    GB_status2("Compressing %i sequences in tree", seqcount);
520                    GB_status(0);
521
522                    for (si=0;si<seqcount;si++){
523                        int        mi     = seqs[si].master;
524                        GB_Master *master = masters[mi];
525                        GBDATA    *gbd    = seqs[si].gbd;
526
527                        if (GB_read_clock(gbd) >= main_clock){
528                            GB_warning("A species seems to be more than once in the tree");
529                        }
530                        else {
531                            char *seq        = GB_read_string(gbd);
532                            int   seq_len    = GB_read_string_count(gbd);
533                            long  sizen      = GB_read_memuse(gbd);
534                            char *seqm       = GB_read_string(master->gbd);
535                            int   master_len = GB_read_string_count(master->gbd);
536                            long  sizes;
537                            char *ss         = gb_compress_sequence_by_master(gbd,seqm,master_len,mi,ali_quark,seq,seq_len,&sizes);
538
539                            gb_write_compressed_pntr(gbd,ss,sizes,seq_len);
540                            sizes = GB_read_memuse(gbd); // check real usage
541
542                            sumnew += sizes;
543                            sumold += sizen;
544                            sumorg += seq_len;
545
546                            free(seqm);
547                            free(seq);
548                        }
549
550                        if (GB_status((si+1)/(double)seqcount)) {
551                            error = "User abort";
552                            break;
553                        }
554                    }
555                }
556
557                /* Compress rest of sequences */
558                if (!error) {
559                    int pass; /* pass 1 : count species to compress, pass 2 : compress species */
560                    int speciesNotInTree = 0;
561
562                    GB_status2("Compressing sequences NOT in tree");
563                    GB_status(0);
564
565                    for (pass = 1; pass <= 2; ++pass) {
566                        GBDATA *gb_species_data = GB_search(gb_main,"species_data",GB_CREATE_CONTAINER);
567                        GBDATA *gb_species;
568                        int     count           = 0;
569
570                        for (gb_species = GBT_first_species_rel_species_data(gb_species_data);
571                             gb_species;
572                             gb_species = GBT_next_species(gb_species))
573                        {
574                            GBDATA *gbd = GBT_read_sequence(gb_species,ali_name);
575
576                            if(!gbd) continue;
577                            if (GB_read_clock(gbd) >= main_clock) continue; /* Compress only those which are not compressed by masters */
578                            count++;
579                            if (pass == 2) {
580                                char *data    = GB_read_string(gbd);
581                                long  seq_len = GB_read_string_count(gbd);
582                                long  size    = GB_read_memuse(gbd);
583
584                                GB_write_string(gbd,"");
585                                GB_write_string(gbd,data);
586                                free(data);
587
588                                sumold += size;
589
590                                size = GB_read_memuse(gbd);
591                                sumnew += size;
592                                sumorg += seq_len;
593
594                                if (GB_status(count/(double)speciesNotInTree)) {
595                                    error = "User abort";
596                                    break;
597                                }
598                            }
599                        }
600                        if (pass == 1) {
601                            speciesNotInTree = count;
602                            if (GB_status2("Compressing %i sequences NOT in tree", speciesNotInTree)) {
603                                error = "User abort";
604                                break;
605                            }
606                        }
607                    }
608                }
609
610                if (!error) {
611                    GB_status2("Compressing %i master-sequences", mastercount);
612                    GB_status(0);
613
614                    /* Compress all masters */
615                    for (si=0;si<mastercount;si++){
616                        int mi = masters[si]->master;
617
618                        if (mi>0) { /*  master available */
619                            GBDATA *gbd = masters[si]->gbd;
620
621                            ad_assert(mi>si); /* we don't want a recursion, because we cannot uncompress sequence compressed masters, Main->gb_master_data is wrong */
622
623                            if (gb_read_nr(gbd) != si) { /* Check database */
624                                GB_internal_error("Sequence Compression: Master Index Conflict");
625                                error = GB_export_error("Sequence Compression: Master Index Conflict");
626                                break;
627                            }
628
629                            {
630                                GB_Master *master     = masters[mi];
631                                char      *seqm       = GB_read_string(master->gbd);
632                                int        master_len = GB_read_string_count(master->gbd);
633                                char      *seq        = GB_read_string(gbd);
634                                int        seq_len    = GB_read_string_count(gbd);
635                                long       sizes;
636                                char      *ss         = gb_compress_sequence_by_master(gbd,seqm,master_len,mi,ali_quark,seq,seq_len,&sizes);
637
638                                gb_write_compressed_pntr(gbd,ss,sizes,seq_len);
639                                sumnew += sizes;
640
641                                free(seq);
642                                free(seqm);
643                            }
644
645                            if (GB_status((si+1)/(double)mastercount)) {
646                                error = "User abort";
647                                break;
648                            }
649                        }
650                        else { // count size of top master
651                            GBDATA *gbd  = masters[si]->gbd;
652                            sumnew      += GB_read_memuse(gbd);
653                        }
654                    }
655
656                    // count size of old master data
657                    if (!error) {
658                        GBDATA *gb_omaster;
659                        for (gb_omaster = GB_entry(old_gb_master_ali, "@master");
660                             gb_omaster;
661                             gb_omaster = GB_nextEntry(gb_omaster))
662                        {
663                            long size  = GB_read_memuse(gb_omaster);
664                            sumold    += size;
665                        }
666                    }
667
668                    if (!error) {
669                        char *sizeOrg = strdup(GBS_readable_size(sumorg));
670                        char *sizeOld = strdup(GBS_readable_size(sumold));
671                        char *sizeNew = strdup(GBS_readable_size(sumnew));
672
673                        GB_warningf("Alignment '%s':\n"
674                                    "    Uncompressed data:   %7s\n"
675                                    "    Old compressed data: %7s = %6.2f%%\n"
676                                    "    New compressed data: %7s = %6.2f%%",
677                                    ali_name, sizeOrg,
678                                    sizeOld, (100.0*sumold)/sumorg,
679                                    sizeNew, (100.0*sumnew)/sumorg);
680
681                        free(sizeNew);
682                        free(sizeOld);
683                        free(sizeOrg);
684                    }
685                }
686
687
688                if (!error) {
689                    if (old_gb_master_ali){
690                        error = GB_delete(old_gb_master_ali);
691                    }
692                    Main->keys[ali_quark].gb_master_ali = gb_master_ali;
693                }
694
695                // free data
696                free(seqs);
697                for (si=0;si<mastercount;si++) free(masters[si]);
698                free(masters);
699            }
700        }
701    }
702    return error;
703}
704
705/* Compress sequences, call only outside a transaction */
706GB_ERROR GBT_compress_sequence_tree2(GBDATA *gb_main, const char *tree_name, const char *ali_name){
707    GB_ERROR      error     = 0;
708    GB_CTREE     *ctree;
709    GB_UNDO_TYPE  undo_type = GB_get_requested_undo_type(gb_main);
710    GB_MAIN_TYPE *Main      = GB_MAIN(gb_main);
711    char         *to_free   = 0;
712
713    if (Main->transaction>0){
714        GB_internal_error("Internal Error: Compress Sequences called during a running transaction");
715        return GB_export_error("Internal Error: Compress Sequences called during a running transaction");
716    }
717
718    GB_request_undo_type(gb_main,GB_UNDO_KILL);
719    GB_begin_transaction(gb_main);
720    GB_push_my_security(gb_main);
721
722    if (!tree_name || !strlen(tree_name)) {
723        to_free   = GBT_find_largest_tree(gb_main);
724        tree_name = to_free;
725    }
726    ctree = (GB_CTREE *)GBT_read_tree(gb_main,(char *)tree_name,-sizeof(GB_CTREE));
727    if (!ctree) {
728        error = GB_export_errorf("Tree %s not found in database",tree_name);
729    }
730    else {
731        error             = GBT_link_tree((GBT_TREE *)ctree, gb_main, GB_FALSE, 0, 0);
732        if (!error) error = compress_sequence_tree(gb_main,ctree,ali_name);
733        GBT_delete_tree((GBT_TREE *)ctree);
734    }
735
736    GB_pop_my_security(gb_main);
737    if (error) {
738        GB_abort_transaction(gb_main);
739    }
740    else {
741        GB_commit_transaction(gb_main);
742        GB_disable_quicksave(gb_main,"Database optimized");
743    }
744    GB_request_undo_type(gb_main,undo_type);
745
746    if (to_free) free(to_free);
747
748#if defined(SAVE_COMPRESSION_TREE_TO_DB)
749    error = "fake error";
750#endif /* SAVE_COMPRESSION_TREE_TO_DB */
751
752    return error;
753}
754
755void GBT_compression_test(void *dummy, GBDATA *gb_main) {
756    GB_ERROR  error     = GB_begin_transaction(gb_main);
757    char     *ali_name  = GBT_get_default_alignment(gb_main);
758    char     *tree_name = GBT_read_string(gb_main, "focus/tree_name");
759
760    GBUSE(dummy);
761    if (!ali_name || !tree_name) error = GB_await_error();
762
763    error = GB_end_transaction(gb_main, error);
764
765    if (!error) {
766        printf("Recompression data in alignment '%s' using tree '%s'\n", ali_name, tree_name);
767        error = GBT_compress_sequence_tree2(gb_main, tree_name, ali_name);
768    }
769
770    if (error) GB_warning(error);
771    free(tree_name);
772    free(ali_name);
773}
774
775
776
777/* ******************** Decompress Sequences ******************** */
778
779char *g_b_uncompress_single_sequence_by_master(const char *s, const char *master, long size, long *new_size) {
780    const signed char *source = (signed char *)s;
781    char              *dest;
782    const char        *m      = master;
783    unsigned int       c;
784    int                j;
785    int                i;
786    char              *buffer;
787
788    dest = buffer = GB_give_other_buffer((char *)source,size);
789
790    for (i=size;i;) {
791        j = *(source++);
792        if (j>0) {      /* uncompressed data block */
793            if (j>i) j=i;
794            i -= j;
795            for (;j;j--) {
796                c = *(source++);
797                if (!c) c = *m;
798                *(dest++) = c;
799                m++;
800            }
801        }else{          /* equal bytes compressed */
802            if (!j) break;  /* end symbol */
803            if (j== -122) {
804                j = *(source++) & 0xff;
805                j |= ((*(source++)) <<8) &0xff00;
806                j = -j;
807            }
808            c = *(source++);
809            i += j;
810            if (i<0) {
811                GB_internal_error("Internal Error: Missing end in data");
812                j += -i;
813                i = 0;
814            }
815            if (c==0) {
816                memcpy(dest, m, -j);
817                dest+= -j;
818                m+= -j;
819            } else {
820                memset(dest, c, -j);
821                dest+= -j;
822                m+= -j;
823            }
824        }
825    }
826    *(dest++) = 0;              /* NULL of NULL terminated string */
827
828    *new_size = dest-buffer;
829    ad_assert(size == *new_size); // buffer overflow
830
831    return buffer;
832}
833
834char *gb_uncompress_by_sequence(GBDATA *gbd, const char *ss,long size, GB_ERROR *error, long *new_size) {
835    char *dest = 0;
836
837    *error = 0;
838    if (!GB_FATHER(gbd)) {
839        *error = "Cannot uncompress this sequence: Sequence has no father";
840    }
841    else {
842        GB_MAIN_TYPE *Main    = GB_MAIN(gbd);
843        GBDATA       *gb_main = (GBDATA*)Main->data;
844        char         *to_free = GB_check_out_buffer(ss); /* Remove 'ss' from memory management, otherwise load_single_key_data() may destroy it */
845        int           index;
846        GBQUARK       quark;
847
848        {
849            const unsigned char *s = (const unsigned char *)ss;
850
851            index = g_b_read_number2(&s);
852            quark = g_b_read_number2(&s);
853
854            ss = (const char *)s;
855        }
856
857        if (!Main->keys[quark].gb_master_ali){
858            gb_load_single_key_data(gb_main,quark);
859        }
860
861        if (!Main->keys[quark].gb_master_ali){
862            *error = "Cannot uncompress this sequence: Cannot find a master sequence";
863        }
864        else {
865            GBDATA *gb_master = gb_find_by_nr(Main->keys[quark].gb_master_ali,index);
866            if (gb_master){
867                const char *master = GB_read_char_pntr(gb_master); /* make sure that this is not a buffer !!! */
868
869                ad_assert((GB_read_string_count(gb_master)+1) == size); // size mismatch between master and slave
870                dest = g_b_uncompress_single_sequence_by_master(ss, master, size, new_size);
871            }
872            else {
873                *error = GB_await_error();
874            }
875        }
876        free(to_free);
877    }
878
879    return dest;
880}
Note: See TracBrowser for help on using the repository browser.