source: tags/ms_r17q3/ARBDB/adseqcompr.cxx

Last change on this file was 15513, checked in by westram, 9 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 40.4 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : adseqcompr.cxx                                    //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include <arb_progress.h>
12#include <arb_file.h>
13#include <arb_misc.h>
14#include <arb_diff.h>
15#include "ad_cb.h"
16#include "gb_key.h"
17#include "TreeNode.h"
18
19#include <climits>
20
21// --------------------------------------------------------------------------------
22
23#define MAX_SEQUENCE_PER_MASTER 50 // was 18 till May 2008
24
25#if defined(DEBUG)
26// don't do optimize, only create tree and save to DB
27// #define SAVE_COMPRESSION_TREE_TO_DB
28#endif // DEBUG
29
30// --------------------------------------------------------------------------------
31
32struct CompressionRoot : public TreeRoot {
33    CompressionRoot();
34    TreeNode *makeNode() const OVERRIDE;
35    void destroyNode(TreeNode *node) const OVERRIDE;
36};
37
38class CompressionTree FINAL_TYPE : public TreeNode {
39protected:
40    ~CompressionTree() OVERRIDE {}
41    friend class CompressionRoot;
42public:
43
44    // members initialized by init_indices_and_count_sons
45    int index; // master(inner nodes) or sequence(leaf nodes) index
46    int sons;  // sons with sequence or masters (in subtree)
47
48    CompressionTree(CompressionRoot *croot) : TreeNode(croot) {}
49
50    unsigned get_leaf_count() const OVERRIDE {
51        gb_assert(0); // @@@ impl (see also GBT_count_leafs and AP_pos_var::getsize)
52        return 0;
53    }
54    void compute_tree() OVERRIDE {}
55
56    DEFINE_TREE_ACCESSORS(CompressionRoot, CompressionTree);
57};
58
59CompressionRoot::CompressionRoot() : TreeRoot(true) {}
60TreeNode *CompressionRoot::makeNode() const { return new CompressionTree(const_cast<CompressionRoot*>(this)); }
61void CompressionRoot::destroyNode(TreeNode *node) const { delete DOWNCAST(CompressionTree*,node); }
62
63struct Consensus {
64    int            len;
65    char           used[256];
66    unsigned char *con[256];
67};
68
69struct Sequence {
70    GBENTRY *gb_seq;
71    int      master;
72};
73
74struct MasterSequence {
75    GBENTRY *gb_mas;
76    int      master;
77};
78
79// --------------------------------------------------------------------------------
80
81static Consensus *g_b_new_Consensus(long len) {
82    Consensus     *gcon = ARB_calloc<Consensus>(1);
83    unsigned char *data = ARB_calloc<unsigned char>(256*len);
84
85    gcon->len = len;
86
87    for (int i=0; i<256; i++) {
88        gcon->con[i] = data + len*i;
89    }
90    return gcon;
91}
92
93
94static void g_b_delete_Consensus(Consensus *gcon) {
95    free(gcon->con[0]);
96    free(gcon);
97}
98
99
100static void g_b_Consensus_add(Consensus *gcon, unsigned char *seq, long seq_len) {
101    const int max_priority = 255/MAX_SEQUENCE_PER_MASTER;      // No overflow possible
102    gb_assert(max_priority >= 1);
103
104    if (seq_len > gcon->len) seq_len = gcon->len;
105
106    // Search for runs
107    unsigned char *s = seq;
108    int last = 0;
109    int i;
110    int li;
111    int c;
112
113    for (li = i = 0; i < seq_len; i++) {
114        c = *(s++);
115        if (c == last) {
116            continue;
117        }
118        else {
119          inc_hits :
120            int eq_count = i-li;
121            gcon->used[c] = 1;
122            unsigned char *p = gcon->con[last];
123            last = c;
124
125            if (eq_count <= GB_RUNLENGTH_SIZE) {
126                c = max_priority;
127                while (li < i) p[li++] += c;
128            }
129            else {
130                c = max_priority * (GB_RUNLENGTH_SIZE) / eq_count;
131                if (c) {
132                    while (li < i) p[li++] += c;
133                }
134                else {
135                    while (li < i) p[li++] |= 1;
136                }
137            }
138        }
139    }
140    if (li<seq_len) {
141        c = last;
142        i = seq_len;
143        goto inc_hits;
144    }
145}
146
147static char *g_b_Consensus_get_sequence(Consensus *gcon) {
148    int pos;
149
150    unsigned char *s;
151
152    unsigned char *max = ARB_calloc<unsigned char>(gcon->len);
153    char          *seq = ARB_calloc<char>(gcon->len+1);
154
155    memset(seq, '@', gcon->len);
156
157    for (int c = 1; c<256; c++) { // Find maximum frequency of non run
158        if (!gcon->used[c]) continue;
159        s = gcon->con[c];
160        for (pos = 0; pos<gcon->len; pos++) {
161            if (*s > max[pos]) {
162                max[pos] = *s;
163                seq[pos] = c;
164            }
165            s++;
166        }
167    }
168    free(max);
169    return seq;
170}
171
172
173static int g_b_count_leafs(CompressionTree *node) {
174    if (node->is_leaf) return 1;
175    node->gb_node = 0;
176    return (g_b_count_leafs(node->get_leftson()) + g_b_count_leafs(node->get_rightson()));
177}
178
179static void g_b_put_sequences_in_container(CompressionTree *ctree, Sequence *seqs, MasterSequence **masters, Consensus *gcon) {
180    if (ctree->is_leaf) {
181        if (ctree->index >= 0) {
182            GB_CSTR data = GB_read_char_pntr(seqs[ctree->index].gb_seq);
183            long    len  = GB_read_string_count(seqs[ctree->index].gb_seq);
184            g_b_Consensus_add(gcon, (unsigned char *)data, len);
185        }
186    }
187    else if (ctree->index<0) {
188        g_b_put_sequences_in_container(ctree->get_leftson(), seqs, masters, gcon);
189        g_b_put_sequences_in_container(ctree->get_rightson(), seqs, masters, gcon);
190    }
191    else {
192        GB_CSTR data = GB_read_char_pntr(masters[ctree->index]->gb_mas);
193        long    len  = GB_read_string_count(masters[ctree->index]->gb_mas);
194        g_b_Consensus_add(gcon, (unsigned char *)data, len);
195    }
196}
197
198static void g_b_create_master(CompressionTree *node, Sequence *seqs, MasterSequence **masters, int my_master, const char *ali_name, long seq_len, arb_progress& progress) {
199    if (node->is_leaf) {
200        if (node->index >= 0) {
201            GBDATA *gb_data = GBT_find_sequence(node->gb_node, ali_name);
202
203            seqs[node->index].gb_seq = gb_data->as_entry();
204            seqs[node->index].master = my_master;
205        }
206    }
207    else {
208        if (progress.aborted()) return;
209       
210        if (node->index>=0) {
211            masters[node->index]->master = my_master;
212            my_master = node->index;
213        }
214        g_b_create_master(node->get_leftson(), seqs, masters, my_master, ali_name, seq_len, progress);
215        g_b_create_master(node->get_rightson(), seqs, masters, my_master, ali_name, seq_len, progress);
216        if (node->index>=0 && !progress.aborted()) { // build me
217            char      *data;
218            Consensus *gcon = g_b_new_Consensus(seq_len);
219
220            g_b_put_sequences_in_container(node->get_leftson(), seqs, masters, gcon);
221            g_b_put_sequences_in_container(node->get_rightson(), seqs, masters, gcon);
222
223            data = g_b_Consensus_get_sequence(gcon);
224
225            GB_write_string(masters[node->index]->gb_mas, data);
226            GB_write_security_write(masters[node->index]->gb_mas, 7);
227
228            g_b_delete_Consensus(gcon);
229            free(data);
230
231            ++progress;
232        }
233    }
234}
235
236// -------------------------------------
237//      distribute master sequences
238
239static void subtract_sons_from_tree(CompressionTree *node, int subtract) {
240    while (node) {
241        node->sons -= subtract;
242        node        = node->get_father();
243    }
244}
245
246static int set_masters_with_sons(CompressionTree *node, int wantedSons, int *mcount) {
247    if (!node->is_leaf) {
248        if (node->sons == wantedSons) {
249            // insert new master
250            gb_assert(node->index == -1);
251            node->index = *mcount;
252            (*mcount)++;
253
254            subtract_sons_from_tree(node->get_father(), node->sons-1);
255            node->sons = 1;
256        }
257        else if (node->sons>wantedSons) {
258            int lMax = set_masters_with_sons(node->get_leftson(),  wantedSons, mcount);
259            int rMax = set_masters_with_sons(node->get_rightson(), wantedSons, mcount);
260
261            int maxSons = lMax<rMax ? rMax : lMax;
262            if (node->sons <= MAX_SEQUENCE_PER_MASTER && node->sons>maxSons) {
263                maxSons = node->sons;
264            }
265            return maxSons;
266        }
267    }
268    return node->sons <= MAX_SEQUENCE_PER_MASTER ? node->sons : 0;
269}
270
271static int maxCompressionSteps(CompressionTree *node) {
272    if (node->is_leaf) {
273        return 0;
274    }
275
276    int left  = maxCompressionSteps(node->get_leftson());
277    int right = maxCompressionSteps(node->get_rightson());
278
279#if defined(SAVE_COMPRESSION_TREE_TO_DB)
280    freenull(node->name);
281    if (node->index2 != -1) {
282        node->name = GBS_global_string_copy("master_%03i", node->index2);
283    }
284#endif // SAVE_COMPRESSION_TREE_TO_DB
285
286    return (left>right ? left : right) + (node->index == -1 ? 0 : 1);
287}
288
289static int init_indices_and_count_sons(CompressionTree *node, int *scount, const char *ali_name) {
290    if (node->is_leaf) {
291        if (node->gb_node == 0 || !GBT_find_sequence(node->gb_node, (char *)ali_name)) {
292            node->index = -1;
293            node->sons  = 0;
294        }
295        else {
296            node->index = *scount;
297            node->sons  = 1;
298            (*scount)++;
299        }
300    }
301    else {
302        node->index = -1;
303        node->sons  =
304            init_indices_and_count_sons(node->get_leftson(), scount, ali_name) +
305            init_indices_and_count_sons(node->get_rightson(), scount, ali_name);
306    }
307    return node->sons;
308}
309
310static void distribute_masters(CompressionTree *tree, int *mcount, int *max_masters) {
311    int wantedSons = MAX_SEQUENCE_PER_MASTER;
312    while (wantedSons >= 2) {
313        int maxSons = set_masters_with_sons(tree, wantedSons, mcount);
314        wantedSons = maxSons;
315    }
316    gb_assert(tree->sons == 1);
317
318    gb_assert(tree->index != -1);
319    *max_masters = maxCompressionSteps(tree);
320}
321
322// --------------------------------------------------------------------------------
323
324#define MAX_NUMBER 0x7fffffff
325STATIC_ASSERT(MAX_NUMBER <= INT_MAX); // ensure 32-bit-version compatibility!
326
327inline int g_b_read_number2(const unsigned char*& s) {
328    int result;
329    unsigned c0 = *s++;
330    if (c0 & 0x80) {
331        unsigned c1 = *s++;
332        if (c0 & 0x40) {
333            unsigned c2 = *s++;
334            if (c0 & 0x20) {
335                unsigned c3 = *s++;
336                if (c0 & 0x10) {
337                    unsigned c4 = *s++;
338                    result = c4 | (c3<<8) | (c2<<16) | (c1<<24);
339                }
340                else {
341                    result = c3 | (c2<<8) | (c1<<16) | ((c0 & 0x0f)<<24);
342                }
343            }
344            else {
345                result = c2 | (c1<<8) | ((c0 & 0x1f)<<16);
346            }
347        }
348        else {
349            result = c1 | ((c0 & 0x3f)<<8);
350        }
351    }
352    else {
353        result = c0;
354    }
355    gb_assert(result >= 0 && result <= MAX_NUMBER);
356    return result;
357}
358
359inline void g_b_put_number2(int i, unsigned char*& s) {
360    gb_assert(i >= 0 && i <= MAX_NUMBER);
361
362    if (i< 0x80) {
363        *s++ = i;
364    }
365    else {
366        int j;
367        if (i<0x4000) {
368            j = (i>>8) | 0x80; *s++ = j;
369            *s++ = i;
370        }
371        else if (i<0x200000) {
372            j = (i>>16) | 0xC0; *s++ = j;
373            j = (i>>8);         *s++ = j;
374            *s++ = i;
375        }
376        else if (i<0x10000000) {
377            j = (i>>24) | 0xE0; *s++ = j;
378            j = (i>>16);        *s++ = j;
379            j = (i>>8);         *s++ = j;
380            *s++ = i;
381        }
382        else {
383            *s++ = 0xF0;
384            j = (i>>24); *s++ = j;
385            j = (i>>16); *s++ = j;
386            j = (i>>8);  *s++ = j;
387            *s++ = i;
388        }
389    }
390}
391
392// --------------------------------------------------------------------------------
393
394#ifdef UNIT_TESTS
395#ifndef TEST_UNIT_H
396#include <test_unit.h>
397#endif
398
399static arb_test::match_expectation put_read_num_using_bytes(int num_written, int bytes_expected, unsigned char *buffer_expected = NULL) {
400    const int     BUFSIZE = 6;
401    unsigned char buffer[BUFSIZE];
402
403    using namespace arb_test;
404
405    unsigned char INIT = 0xaa;
406    memset(buffer, INIT, BUFSIZE);
407
408    expectation_group expected;
409
410    {
411        unsigned char *bp = buffer;
412        g_b_put_number2(num_written, bp);
413
414        size_t bytes_written = bp-buffer;
415        expected.add(that(bytes_written).is_equal_to(bytes_expected));
416
417        if (buffer_expected) {
418            expected.add(that(arb_test::memory_is_equal(buffer, buffer_expected, bytes_expected)).is_equal_to(true));
419        }
420    }
421    {
422        const unsigned char *bp = buffer;
423
424        int num_read = g_b_read_number2(bp);
425        expected.add(that(num_read).is_equal_to(num_written));
426
427        size_t bytes_read = bp-buffer;
428        expected.add(that(bytes_read).is_equal_to(bytes_expected));
429    }
430
431    expected.add(that(buffer[bytes_expected]).is_equal_to(INIT)); 
432
433    return all().ofgroup(expected);
434}
435
436#define TEST_PUT_READ_NUMBER(num,expect_bytes)         TEST_EXPECTATION(put_read_num_using_bytes(num, expect_bytes))
437#define TEST_PUT_READ_NUMBER__BROKEN(num,expect_bytes) TEST_EXPECTATION__BROKEN(put_read_num_using_bytes(num, expect_bytes))
438
439#define TEST_PUT_NUMBER_BINARY1(num, byte1) do {                \
440        unsigned char buf[1];                                   \
441        buf[0] = byte1;                                         \
442        TEST_EXPECTATION(put_read_num_using_bytes(num, 1, buf));     \
443    } while(0)
444
445#define TEST_PUT_NUMBER_BINARY2(num, byte1, byte2) do {         \
446        unsigned char buf[2];                                   \
447        buf[0] = byte1;                                         \
448        buf[1] = byte2;                                         \
449        TEST_EXPECTATION(put_read_num_using_bytes(num, 2, buf));     \
450    } while(0)
451
452#define TEST_PUT_NUMBER_BINARY3(num, byte1, byte2, byte3) do {  \
453        unsigned char buf[3];                                   \
454        buf[0] = byte1;                                         \
455        buf[1] = byte2;                                         \
456        buf[2] = byte3;                                         \
457        TEST_EXPECTATION(put_read_num_using_bytes(num, 3, buf));     \
458    } while(0)
459
460#define TEST_PUT_NUMBER_BINARY4(num, byte1, byte2, byte3, byte4) do {   \
461        unsigned char buf[4];                                           \
462        buf[0] = byte1;                                                 \
463        buf[1] = byte2;                                                 \
464        buf[2] = byte3;                                                 \
465        buf[3] = byte4;                                                 \
466        TEST_EXPECTATION(put_read_num_using_bytes(num, 4, buf));             \
467    } while(0)
468
469#define TEST_PUT_NUMBER_BINARY5(num, byte1, byte2, byte3, byte4, byte5) do { \
470        unsigned char buf[5];                                           \
471        buf[0] = byte1;                                                 \
472        buf[1] = byte2;                                                 \
473        buf[2] = byte3;                                                 \
474        buf[3] = byte4;                                                 \
475        buf[4] = byte5;                                                 \
476        TEST_EXPECTATION(put_read_num_using_bytes(num, 5, buf));             \
477    } while(0)
478   
479void TEST_put_read_number() {
480    // test that put and read are compatible:
481    TEST_PUT_READ_NUMBER(0x0, 1);
482
483    TEST_PUT_READ_NUMBER(0x7f, 1);
484    TEST_PUT_READ_NUMBER(0x80, 2);
485
486    TEST_PUT_READ_NUMBER(0x3fff, 2);
487    TEST_PUT_READ_NUMBER(0x4000, 3);
488
489    TEST_PUT_READ_NUMBER(0x1fffff, 3);
490    TEST_PUT_READ_NUMBER(0x200000, 4);
491
492    TEST_PUT_READ_NUMBER(0xfffffff, 4);
493
494    TEST_PUT_READ_NUMBER(0x10000000, 5);
495    TEST_PUT_READ_NUMBER(0x7fffffff, 5);
496
497    // test binary compatibility:
498    // (code affects DB content, cannot be changed)
499   
500    TEST_PUT_NUMBER_BINARY1(0x0,  0x00);
501    TEST_PUT_NUMBER_BINARY1(0x7f, 0x7f);
502
503    TEST_PUT_NUMBER_BINARY2(0x80,   0x80, 0x80);
504    TEST_PUT_NUMBER_BINARY2(0x81,   0x80, 0x81);
505    TEST_PUT_NUMBER_BINARY2(0x3fff, 0xbf, 0xff);
506
507    TEST_PUT_NUMBER_BINARY3(0x4000,   0xc0, 0x40, 0x00);
508    TEST_PUT_NUMBER_BINARY3(0x1fffff, 0xdf, 0xff, 0xff);
509
510    TEST_PUT_NUMBER_BINARY4(0x200000,  0xe0, 0x20, 0x00, 0x00);
511    TEST_PUT_NUMBER_BINARY4(0xfffffff, 0xef, 0xff, 0xff, 0xff);
512   
513    TEST_PUT_NUMBER_BINARY5(0x10000000, 0xf0, 0x10, 0x00, 0x00, 0x00);
514    TEST_PUT_NUMBER_BINARY5(0x7fffffff, 0xf0, 0x7f, 0xff, 0xff, 0xff);
515}
516
517#endif // UNIT_TESTS
518
519// --------------------------------------------------------------------------------
520
521
522static char *gb_compress_seq_by_master(const char *master, size_t master_len, int master_index,
523                                       GBQUARK q, const char *seq, size_t seq_len,
524                                       size_t *memsize, int old_flag) {
525    unsigned char *buffer;
526    int            rest = 0;
527    unsigned char *d;
528    int            i, cs, cm;
529    int            last;
530    long           len  = seq_len;
531
532    d = buffer = (unsigned char *)GB_give_other_buffer(seq, seq_len);
533
534    if (seq_len > master_len) {
535        rest = seq_len - master_len;
536        len = master_len;
537    }
538
539    last = -1000;       // Convert Sequence relative to Master
540    for (i = len; i>0; i--) {
541        cm = *(master++);
542        cs = *(seq++);
543        if (cm==cs && cs != last) {
544            *(d++) = 0;
545            last = 1000;
546        }
547        else {
548            *(d++) = cs;
549            last = cs;
550        }
551    }
552    for (i = rest; i>0; i--) {
553        *(d++) = *(seq++);
554    }
555
556    {               // Append run length compression method
557        unsigned char *buffer2;
558        unsigned char *dest2;
559        buffer2 = dest2 = (unsigned char *)GB_give_other_buffer((char *)buffer, seq_len+100);
560        *(dest2++) = GB_COMPRESSION_SEQUENCE | old_flag;
561
562        g_b_put_number2(master_index, dest2); // Tags
563        g_b_put_number2(q, dest2);
564
565        gb_compress_equal_bytes_2((char *)buffer, seq_len, memsize, (char *)dest2); // append runlength compressed sequences to tags
566
567        *memsize = *memsize + (dest2-buffer2);
568        return (char *)buffer2;
569    }
570}
571
572static char *gb_compress_sequence_by_master(GBDATA *gbd, const char *master, size_t master_len, int master_index,
573                                            GBQUARK q, const char *seq, size_t seq_len, size_t *memsize)
574{
575    size_t  size;
576    char   *is  = gb_compress_seq_by_master(master, master_len, master_index, q, seq, seq_len, &size, GB_COMPRESSION_LAST);
577    char   *res = gb_compress_data(gbd, 0, is, size, memsize, ~(GB_COMPRESSION_DICTIONARY|GB_COMPRESSION_SORTBYTES|GB_COMPRESSION_RUNLENGTH), true);
578    return res;
579}
580
581static GB_ERROR compress_sequence_tree(GBCONTAINER *gb_main, CompressionTree *tree, const char *ali_name) {
582    GB_ERROR error      = 0;
583    long     ali_len    = GBT_get_alignment_len(gb_main, ali_name);
584    int      main_clock = GB_read_clock(gb_main);
585
586    GB_ERROR warning = NULL;
587
588    if (ali_len<0) {
589        warning = GBS_global_string("Skipping alignment '%s' (not a valid alignment; len=%li).", ali_name, ali_len);
590        GB_clear_error();
591    }
592    else {
593        int leafcount = g_b_count_leafs(tree);
594        if (!leafcount) {
595            error = "Tree is empty";
596        }
597        else {
598            arb_progress tree_progress("Compressing sequences", 4);
599           
600            // Distribute masters in tree
601            int mastercount   = 0;
602            int max_compSteps = 0; // in one branch
603            int seqcount      = 0;
604
605            init_indices_and_count_sons(tree, &seqcount, ali_name);
606            if (!seqcount) {
607                warning = GBS_global_string("Tree contains no sequences with data in '%s'\n"
608                                            "Skipping compression for this alignment",
609                                            ali_name);
610            }
611            else {
612                distribute_masters(tree, &mastercount, &max_compSteps);
613
614#if defined(SAVE_COMPRESSION_TREE_TO_DB)
615                {
616                    error = GBT_link_tree(tree, gb_main, 0, NULL, NULL);
617                    if (!error) error = GBT_write_tree(gb_main, 0, "tree_compression_new", tree);
618                    GB_information("Only generated compression tree (do NOT save DB anymore)");
619                    return error;
620                }
621#endif // SAVE_COMPRESSION_TREE_TO_DB
622
623                // detect degenerated trees
624                {
625                    int min_masters   = ((seqcount-1)/MAX_SEQUENCE_PER_MASTER)+1;
626                    int min_compSteps = 1;
627                    {
628                        int m = min_masters;
629                        while (m>1) {
630                            m            = ((m-1)/MAX_SEQUENCE_PER_MASTER)+1;
631                            min_masters += m;
632                            min_compSteps++;
633                        }
634                    }
635
636                    int acceptable_masters   = (3*min_masters)/2; // accept 50% overhead
637                    int acceptable_compSteps = 11*min_compSteps; // accept 1000% overhead
638
639                    if (mastercount>acceptable_masters || max_compSteps>acceptable_compSteps) {
640                        GB_warningf("Tree is ill-suited for compression (cause of deep branches)\n"
641                                    "                    Used tree   Optimal tree   Overhead\n"
642                                    "Compression steps       %5i          %5i      %4i%% (speed)\n"
643                                    "Master sequences        %5i          %5i      %4i%% (size)\n"
644                                    "If you like to restart with a better tree,\n"
645                                    "press 'Abort' to stop compression",
646                                    max_compSteps, min_compSteps, (100*max_compSteps)/min_compSteps-100,
647                                    mastercount, min_masters, (100*mastercount)/min_masters-100);
648                    }
649                }
650
651                gb_assert(mastercount>0);
652            }
653
654            if (!warning) {
655                GBCONTAINER        *gb_master_ali     = 0;
656                GBDATA             *old_gb_master_ali = 0;
657                Sequence           *seqs              = 0;
658                GB_MAIN_TYPE       *Main              = GB_MAIN(gb_main);
659                GBQUARK             ali_quark         = gb_find_or_create_quark(Main, ali_name);
660                unsigned long long  sumorg            = 0;
661                unsigned long long  sumold            = 0;
662                unsigned long long  sumnew            = 0;
663                int                 si;
664
665                MasterSequence **masters; ARB_calloc(masters, leafcount);
666
667                {
668                    char *masterfoldername = GBS_global_string_copy("%s/@master_data/@%s", GB_SYSTEM_FOLDER, ali_name);
669                    old_gb_master_ali      = GBDATA::as_container(GB_search(gb_main, masterfoldername, GB_FIND));
670                    free(masterfoldername);
671                }
672
673                // create masters
674                if (!error) {
675                    {
676                        char *master_data_name = GBS_global_string_copy("%s/@master_data", GB_SYSTEM_FOLDER);
677                        char *master_name      = GBS_global_string_copy("@%s", ali_name);
678
679                        GBCONTAINER *gb_master_data = gb_search(gb_main, master_data_name, GB_CREATE_CONTAINER, 1)->as_container();
680
681                        // create a master container, the old is deleted as soon as all sequences are compressed by the new method
682                        gb_master_ali = gb_create_container(gb_master_data, master_name);
683                        GB_write_security_delete(gb_master_ali, 7);
684
685                        free(master_name);
686                        free(master_data_name);
687                    }
688                    for (si = 0; si<mastercount; si++) {
689                        ARB_calloc(masters[si], 1);
690                        masters[si]->gb_mas = gb_create(gb_master_ali, "@master", GB_STRING);
691                    }
692                    ARB_calloc(seqs, leafcount);
693
694                    if (!error) {
695                        arb_progress progress("Building master sequences", mastercount);
696                        g_b_create_master(tree, seqs, masters, -1, ali_name, ali_len, progress);
697                       
698                        error = progress.error_if_aborted();
699                    }
700                }
701                tree_progress.inc_and_check_user_abort(error);
702
703                // Compress sequences in tree
704                if (!error) {
705                    arb_progress progress("Compressing sequences in tree", seqcount);
706
707                    for (si=0; si<seqcount && !error; si++) {
708                        int             mi     = seqs[si].master;
709                        MasterSequence *master = masters[mi];
710                        GBDATA         *gbd    = seqs[si].gb_seq;
711
712                        if (GB_read_clock(gbd) >= main_clock) {
713                            GB_warning("A species seems to be more than once in the tree");
714                        }
715                        else {
716                            char   *seq        = GB_read_string(gbd);
717                            int     seq_len    = GB_read_string_count(gbd);
718                            long    sizen      = GB_read_memuse(gbd);
719                            char   *seqm       = GB_read_string(master->gb_mas);
720                            int     master_len = GB_read_string_count(master->gb_mas);
721                            size_t  sizes;
722                            char   *ss         = gb_compress_sequence_by_master(gbd, seqm, master_len, mi, ali_quark, seq, seq_len, &sizes);
723
724                            gb_write_compressed_pntr(gbd->as_entry(), ss, sizes, seq_len);
725                            sizes = GB_read_memuse(gbd); // check real usage
726
727                            sumnew += sizes;
728                            sumold += sizen;
729                            sumorg += seq_len;
730
731                            free(seqm);
732                            free(seq);
733                        }
734
735                        progress.inc_and_check_user_abort(error);
736                    }
737                }
738                tree_progress.inc_and_check_user_abort(error);
739
740                // Compress rest of sequences
741                if (!error) {
742                    int pass; // pass 1 : count species to compress, pass 2 : compress species
743                    int speciesNotInTree = 0;
744
745                    SmartPtr<arb_progress> progress;
746
747                    for (pass = 1; pass <= 2; ++pass) {
748                        GBDATA *gb_species_data = GBT_get_species_data(gb_main);
749                        GBDATA *gb_species;
750                        int     count           = 0;
751
752                        for (gb_species = GBT_first_species_rel_species_data(gb_species_data);
753                             gb_species;
754                             gb_species = GBT_next_species(gb_species))
755                        {
756                            GBDATA *gbd = GBT_find_sequence(gb_species, ali_name);
757
758                            if (!gbd) continue;
759                            if (GB_read_clock(gbd) >= main_clock) continue; // Compress only those which are not compressed by masters
760                            count++;
761                            if (pass == 2) {
762                                char *data    = GB_read_string(gbd);
763                                long  seq_len = GB_read_string_count(gbd);
764                                long  size    = GB_read_memuse(gbd);
765
766                                GB_write_string(gbd, ""); // force recompression
767                                GB_write_string(gbd, data);
768                                free(data);
769
770                                sumold += size;
771
772                                size = GB_read_memuse(gbd);
773                                sumnew += size;
774                                sumorg += seq_len;
775
776                                progress->inc_and_check_user_abort(error);
777                            }
778                        }
779                        if (pass == 1) {
780                            speciesNotInTree = count;
781                            if (speciesNotInTree>0) {
782                                progress = new arb_progress("Compressing sequences NOT in tree", speciesNotInTree);
783                            }
784                        }
785                    }
786                }
787                tree_progress.inc_and_check_user_abort(error);
788
789                if (!error) {
790                    arb_progress progress("Compressing master-sequences", mastercount);
791
792                    // Compress all masters
793                    for (si=0; si<mastercount; si++) {
794                        int mi = masters[si]->master;
795
796                        if (mi>0) { //  master available
797                            GBDATA *gbd = masters[si]->gb_mas;
798
799                            gb_assert(mi>si); // we don't want a recursion, because we cannot uncompress sequence compressed masters, Main->gb_master_data is wrong
800
801                            if (gb_read_nr(gbd) != si) { // Check database
802                                GB_internal_error("Sequence Compression: Master Index Conflict");
803                                error = GB_export_error("Sequence Compression: Master Index Conflict");
804                                break;
805                            }
806
807                            {
808                                MasterSequence *master     = masters[mi];
809                                char           *seqm       = GB_read_string(master->gb_mas);
810                                int             master_len = GB_read_string_count(master->gb_mas);
811                                char           *seq        = GB_read_string(gbd);
812                                int             seq_len    = GB_read_string_count(gbd);
813                                size_t          sizes;
814                                char           *ss         = gb_compress_sequence_by_master(gbd, seqm, master_len, mi, ali_quark, seq, seq_len, &sizes);
815
816                                gb_write_compressed_pntr(gbd->as_entry(), ss, sizes, seq_len);
817                                sumnew += sizes;
818
819                                free(seq);
820                                free(seqm);
821                            }
822
823                            progress.inc_and_check_user_abort(error);
824                        }
825                        else { // count size of top master
826                            GBDATA *gbd  = masters[si]->gb_mas;
827                            sumnew      += GB_read_memuse(gbd);
828
829                            progress.inc_and_check_user_abort(error);
830                        }
831                    }
832
833                    // count size of old master data
834                    if (!error) {
835                        GBDATA *gb_omaster;
836                        for (gb_omaster = GB_entry(old_gb_master_ali, "@master");
837                             gb_omaster;
838                             gb_omaster = GB_nextEntry(gb_omaster))
839                        {
840                            long size  = GB_read_memuse(gb_omaster);
841                            sumold    += size;
842                        }
843                    }
844
845                    if (!error) {
846                        char *sizeOrg = ARB_strdup(GBS_readable_size(sumorg, "b"));
847                        char *sizeOld = ARB_strdup(GBS_readable_size(sumold, "b"));
848                        char *sizeNew = ARB_strdup(GBS_readable_size(sumnew, "b"));
849
850                        GB_warningf("Alignment '%s':\n"
851                                    "    Uncompressed data:   %7s\n"
852                                    "    Old compressed data: %7s = %6.2f%%\n"
853                                    "    New compressed data: %7s = %6.2f%%",
854                                    ali_name, sizeOrg,
855                                    sizeOld, (100.0*sumold)/sumorg,
856                                    sizeNew, (100.0*sumnew)/sumorg);
857
858                        free(sizeNew);
859                        free(sizeOld);
860                        free(sizeOrg);
861                    }
862                }
863                tree_progress.inc_and_check_user_abort(error);
864
865                if (!error) {
866                    if (old_gb_master_ali) error = GB_delete(old_gb_master_ali);
867                    Main->keys[ali_quark].gb_master_ali = gb_master_ali;
868                }
869
870                // free data
871                free(seqs);
872                for (si=0; si<mastercount; si++) free(masters[si]);
873                free(masters);
874            }
875            else {
876                tree_progress.done();
877            }
878        }
879    }
880
881    if (warning) GB_information(warning);
882
883    return error;
884}
885
886GB_ERROR GBT_compress_sequence_tree2(GBDATA *gbd, const char *tree_name, const char *ali_name) { // goes to header: __ATTR__USERESULT // @@@ rename function
887    // Compress sequences, call only outside a transaction
888    GB_ERROR      error = NULL;
889    GB_MAIN_TYPE *Main  = GB_MAIN(gbd);
890
891    if (Main->get_transaction_level() > 0) {
892        error = "Compress Sequences called while transaction running";
893        GB_internal_error(error);
894    }
895    else {
896        GBCONTAINER  *gb_main        = Main->root_container;
897        GB_UNDO_TYPE  prev_undo_type = GB_get_requested_undo_type(gb_main);
898
899        error = GB_request_undo_type(gb_main, GB_UNDO_KILL);
900        if (!error) {
901            error = GB_begin_transaction(gb_main);
902            if (!error) {
903                GB_push_my_security(gb_main);
904
905                if (!tree_name || !strlen(tree_name)) {
906                    tree_name = GBT_name_of_largest_tree(gb_main);
907                }
908
909                {
910                    CompressionTree *ctree = DOWNCAST(CompressionTree*, GBT_read_tree(gb_main, tree_name, new CompressionRoot));
911                    if (!ctree) error      = GB_await_error();
912                    else {
913                        error             = GBT_link_tree(ctree, gb_main, false, 0, 0);
914                        if (!error) error = compress_sequence_tree(gb_main, ctree, ali_name);
915                        destroy(ctree);
916                    }
917                }
918                if (!error) GB_disable_quicksave(gb_main, "Database optimized");
919
920                GB_pop_my_security(gb_main);
921                error = GB_end_transaction(gb_main, error);
922            }
923            ASSERT_NO_ERROR(GB_request_undo_type(gb_main, prev_undo_type));
924        }
925
926#if defined(SAVE_COMPRESSION_TREE_TO_DB)
927        error = "fake error";
928#endif // SAVE_COMPRESSION_TREE_TO_DB
929    }
930    return error;
931}
932
933#ifdef DEBUG
934
935void GBT_compression_test(struct Unfixed_cb_parameter *, GBDATA *gb_main) {
936    GB_ERROR  error     = GB_begin_transaction(gb_main);
937    char     *ali_name  = GBT_get_default_alignment(gb_main);
938    char     *tree_name = GBT_read_string(gb_main, "focus/tree_name");
939
940    // GBUSE(dummy);
941    if (!ali_name || !tree_name) error = GB_await_error();
942
943    error = GB_end_transaction(gb_main, error);
944
945    if (!error) {
946        printf("Recompression data in alignment '%s' using tree '%s'\n", ali_name, tree_name);
947        error = GBT_compress_sequence_tree2(gb_main, tree_name, ali_name);
948    }
949
950    if (error) GB_warning(error);
951    free(tree_name);
952    free(ali_name);
953}
954
955#endif
956
957// ******************** Decompress Sequences ********************
958
959static char *g_b_uncompress_single_sequence_by_master(const char *s, const char *master, size_t size, size_t *new_size) {
960    const signed char *source = (signed char *)s;
961    char              *dest;
962    const char        *m      = master;
963    unsigned int       c;
964    int                j;
965    int                i;
966    char              *buffer;
967
968    dest = buffer = GB_give_other_buffer((char *)source, size);
969
970    for (i=size; i;) {
971        j = *(source++);
972        if (j>0) {                                  // uncompressed data block
973            if (j>i) j=i;
974            i -= j;
975            for (; j; j--) {
976                c = *(source++);
977                if (!c) c = *m;
978                *(dest++) = c;
979                m++;
980            }
981        }
982        else {                                      // equal bytes compressed
983            if (!j) break;                          // end symbol
984            if (j == -122) {
985                j = *(source++) & 0xff;
986                j |= ((*(source++)) << 8) &0xff00;
987                j = -j;
988            }
989            c = *(source++);
990            i += j;
991            if (i<0) {
992                GB_internal_error("Internal Error: Missing end in data");
993                j += -i;
994                i = 0;
995            }
996            if (c==0) {
997                memcpy(dest, m, -j);
998                dest += -j;
999                m += -j;
1000            }
1001            else {
1002                memset(dest, c, -j);
1003                dest += -j;
1004                m += -j;
1005            }
1006        }
1007    }
1008    *(dest++) = 0;              // NULL of NULL terminated string
1009
1010    *new_size = dest-buffer;
1011    gb_assert(size == *new_size); // buffer overflow
1012
1013    return buffer;
1014}
1015
1016char *gb_uncompress_by_sequence(GBDATA *gbd, const char *ss, size_t size, GB_ERROR *error, size_t *new_size) {
1017    char *dest = 0;
1018
1019    *error = 0;
1020
1021    GB_MAIN_TYPE *Main = gb_get_main_during_cb();
1022    if (!Main && GB_FATHER(gbd)) Main = GB_MAIN(gbd);
1023
1024    if (!Main) {
1025        *error = "Can not uncompress this sequence (neighter has father nor inside callback)";
1026    }
1027    else {
1028        GBDATA  *gb_main = Main->gb_main();
1029        char    *to_free = GB_check_out_buffer(ss);    // Remove 'ss' from memory management, otherwise load_single_key_data() may destroy it
1030        int      index;
1031        GBQUARK  quark;
1032
1033        {
1034            const unsigned char *s = (const unsigned char *)ss;
1035
1036            index = g_b_read_number2(s);
1037            quark = g_b_read_number2(s);
1038
1039            ss = (const char *)s;
1040        }
1041
1042        if (!Main->keys[quark].gb_master_ali) {
1043            gb_load_single_key_data(gb_main, quark);
1044        }
1045
1046        if (!Main->keys[quark].gb_master_ali) {
1047            *error = "Cannot uncompress this sequence: Cannot find a master sequence";
1048        }
1049        else {
1050            GBDATA *gb_master = gb_find_by_nr(Main->keys[quark].gb_master_ali, index);
1051            if (gb_master) {
1052                const char *master = GB_read_char_pntr(gb_master); // make sure that this is not a buffer !!!
1053
1054                gb_assert((GB_read_string_count(gb_master)+1) == size); // size mismatch between master and slave
1055                dest = g_b_uncompress_single_sequence_by_master(ss, master, size, new_size);
1056            }
1057            else {
1058                *error = GB_await_error();
1059            }
1060        }
1061        free(to_free);
1062    }
1063
1064    return dest;
1065}
1066
1067// --------------------------------------------------------------------------------
1068
1069#ifdef UNIT_TESTS
1070#ifndef TEST_UNIT_H
1071#include <test_unit.h>
1072#endif
1073
1074// #define TEST_AUTO_UPDATE // uncomment to auto-update expected result DB
1075
1076void TEST_SLOW_sequence_compression() {
1077    const char *source     = "TEST_nuc.arb";
1078    const char *compressed = "TEST_nuc_seqcompr.arb";
1079    const char *expected   = "TEST_nuc_seqcompr_exp.arb";
1080    const char *aliname    = "ali_16s";
1081
1082    GB_shell shell;
1083
1084    const int  SEQ2COMPARE = 7;
1085    char      *seq_exp[SEQ2COMPARE];
1086
1087    {
1088        GBDATA *gb_main;
1089        TEST_EXPECT_RESULT__NOERROREXPORTED(gb_main = GB_open(source, "rw"));
1090
1091        {
1092            GB_transaction ta(gb_main);
1093            int            count = 0;
1094
1095            for (GBDATA *gb_species = GBT_first_species(gb_main);
1096                 gb_species && count<SEQ2COMPARE;
1097                 gb_species = GBT_next_species(gb_species), ++count)
1098            {
1099                GBDATA *gb_seq = GBT_find_sequence(gb_species, aliname);
1100                seq_exp[count] = GB_read_string(gb_seq);
1101            }
1102        }
1103
1104        TEST_EXPECT_NO_ERROR(GBT_compress_sequence_tree2(gb_main, "tree_nuc", aliname));
1105        TEST_EXPECT_NO_ERROR(GB_save_as(gb_main, compressed, "b"));
1106        GB_close(gb_main);
1107    }
1108#if defined(TEST_AUTO_UPDATE)
1109    TEST_COPY_FILE(compressed, expected);
1110#endif
1111    TEST_EXPECT_FILES_EQUAL(compressed, expected);
1112
1113    {
1114        GBDATA *gb_main;
1115        TEST_EXPECT_RESULT__NOERROREXPORTED(gb_main = GB_open(compressed, "rw"));
1116        {
1117            GB_transaction ta(gb_main);
1118            int            count = 0;
1119
1120            for (GBDATA *gb_species = GBT_first_species(gb_main);
1121                 gb_species && count<SEQ2COMPARE;
1122                 gb_species = GBT_next_species(gb_species), ++count)
1123            {
1124                GBDATA *gb_seq = GBT_find_sequence(gb_species, aliname);
1125                char   *seq    = GB_read_string(gb_seq);
1126
1127                TEST_EXPECT_EQUAL(seq, seq_exp[count]);
1128
1129                freenull(seq_exp[count]);
1130                free(seq);
1131            }
1132        }
1133        GB_close(gb_main);
1134    }
1135
1136    TEST_EXPECT_ZERO_OR_SHOW_ERRNO(GB_unlink(compressed));
1137}
1138TEST_PUBLISH(TEST_SLOW_sequence_compression);
1139
1140#endif // UNIT_TESTS
1141
1142// --------------------------------------------------------------------------------
1143
1144
Note: See TracBrowser for help on using the repository browser.