source: branches/profile/EDIT4/ED4_mini_classes.cxx

Last change on this file was 12238, checked in by westram, 5 years ago
  • use typed CBs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.6 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : ED4_mini_classes.cxx                              //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11
12#include "ed4_class.hxx"
13#include "ed4_edit_string.hxx"
14#include "ed4_awars.hxx"
15#include "ed4_tools.hxx"
16
17#include <aw_awar.hxx>
18#include <aw_root.hxx>
19#include <iupac.h>
20
21#include <cctype>
22
23// ------------------------
24//      ED4_bases_table
25
26#ifdef DEBUG
27// # define COUNT_BASES_TABLE_SIZE
28#endif
29
30#ifdef COUNT_BASES_TABLE_SIZE
31static long bases_table_size = 0L;
32#endif
33
34void ED4_bases_table::init(int length)
35{
36    if (length) {
37        if (no_of_bases.shortTable) {
38            delete no_of_bases.shortTable;
39            no_of_bases.shortTable = 0;
40#ifdef COUNT_BASES_TABLE_SIZE
41            bases_table_size -= (no_of_entries+1)*table_entry_size;
42#endif
43        }
44
45        no_of_entries = length;
46        int new_size = (no_of_entries+1)*table_entry_size;
47        no_of_bases.shortTable = (unsigned char*)new char[new_size];
48        memset(no_of_bases.shortTable, 0, new_size);
49#ifdef COUNT_BASES_TABLE_SIZE
50        bases_table_size += new_size;
51        printf("bases_table_size = %li\n", bases_table_size);
52#endif
53    }
54}
55
56void ED4_bases_table::expand_table_entry_size() { // converts short table into long table
57    e4_assert(table_entry_size==SHORT_TABLE_ELEM_SIZE);
58
59    int *new_table = new int[no_of_entries+1];
60    int i;
61
62    for (i=0; i<=no_of_entries; i++) {
63        new_table[i] = no_of_bases.shortTable[i];
64    }
65    delete [] no_of_bases.shortTable;
66
67#ifdef COUNT_BASES_TABLE_SIZE
68    bases_table_size += (no_of_entries+1)*(LONG_TABLE_ELEM_SIZE-SHORT_TABLE_ELEM_SIZE);
69    printf("bases_table_size = %li\n", bases_table_size);
70#endif
71
72    no_of_bases.longTable = new_table;
73    table_entry_size = LONG_TABLE_ELEM_SIZE;
74}
75
76void ED4_bases_table::add(const ED4_bases_table& other, int start, int end)
77{
78    e4_assert(no_of_entries==other.no_of_entries);
79    e4_assert(start>=0);
80    e4_assert(end<no_of_entries);
81    e4_assert(start<=end);
82
83    int i;
84    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
85        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
86            for (i=start; i<=end; i++) {
87                set_elem_short(i, get_elem_short(i)+other.get_elem_short(i));
88            }
89        }
90        else {
91            for (i=start; i<=end; i++) {
92                set_elem_short(i, get_elem_short(i)+other.get_elem_long(i));
93            }
94        }
95    }
96    else {
97        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
98            for (i=start; i<=end; i++) {
99                set_elem_long(i, get_elem_long(i)+other.get_elem_short(i));
100            }
101        }
102        else {
103            for (i=start; i<=end; i++) {
104                set_elem_long(i, get_elem_long(i)+other.get_elem_long(i));
105            }
106        }
107    }
108}
109void ED4_bases_table::sub(const ED4_bases_table& other, int start, int end)
110{
111    e4_assert(no_of_entries==other.no_of_entries);
112    e4_assert(start>=0);
113    e4_assert(end<no_of_entries);
114    e4_assert(start<=end);
115
116    int i;
117    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
118        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
119            for (i=start; i<=end; i++) {
120                set_elem_short(i, get_elem_short(i)-other.get_elem_short(i));
121            }
122        }
123        else {
124            for (i=start; i<=end; i++) {
125                set_elem_short(i, get_elem_short(i)-other.get_elem_long(i));
126            }
127        }
128    }
129    else {
130        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
131            for (i=start; i<=end; i++) {
132                set_elem_long(i, get_elem_long(i)-other.get_elem_short(i));
133            }
134        }
135        else {
136            for (i=start; i<=end; i++) {
137                set_elem_long(i, get_elem_long(i)-other.get_elem_long(i));
138            }
139        }
140    }
141}
142void ED4_bases_table::sub_and_add(const ED4_bases_table& Sub, const ED4_bases_table& Add, PosRange range)
143{
144    e4_assert(no_of_entries==Sub.no_of_entries);
145    e4_assert(no_of_entries==Add.no_of_entries);
146
147    int start = range.start();
148    int end   = range.end();
149
150    e4_assert(start>=0);
151    e4_assert(end<no_of_entries);
152    e4_assert(start<=end);
153
154    int i;
155    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
156        if (Sub.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
157            if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
158                for (i=start; i<=end; i++) {
159                    set_elem_short(i, get_elem_short(i)-Sub.get_elem_short(i)+Add.get_elem_short(i));
160                }
161            }
162            else {
163                for (i=start; i<=end; i++) {
164                    set_elem_short(i, get_elem_short(i)-Sub.get_elem_short(i)+Add.get_elem_long(i));
165                }
166            }
167        }
168        else {
169            if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
170                for (i=start; i<=end; i++) {
171                    set_elem_short(i, get_elem_short(i)-Sub.get_elem_long(i)+Add.get_elem_short(i));
172                }
173            }
174            else {
175                for (i=start; i<=end; i++) {
176                    set_elem_short(i, get_elem_short(i)-Sub.get_elem_long(i)+Add.get_elem_long(i));
177                }
178            }
179        }
180    }
181    else {
182        if (Sub.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
183            if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
184                for (i=start; i<=end; i++) {
185                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_short(i));
186                }
187            }
188            else {
189                for (i=start; i<=end; i++) {
190                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_long(i));
191                }
192            }
193        }
194        else {
195            if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
196                for (i=start; i<=end; i++) {
197                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_short(i));
198                }
199            }
200            else {
201                for (i=start; i<=end; i++) {
202                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_long(i));
203                }
204            }
205        }
206    }
207}
208
209int ED4_bases_table::firstDifference(const ED4_bases_table& other, int start, int end, int *firstDifferentPos) const
210{
211    int i;
212    int result = 0;
213
214    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
215        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
216            for (i=start; i<=end; i++) {
217                if (get_elem_short(i)!=other.get_elem_short(i)) {
218                    *firstDifferentPos = i;
219                    result = 1;
220                    break;
221                }
222            }
223        }
224        else {
225            for (i=start; i<=end; i++) {
226                if (get_elem_short(i)!=other.get_elem_long(i)) {
227                    *firstDifferentPos = i;
228                    result = 1;
229                    break;
230                }
231            }
232        }
233    }
234    else {
235        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
236            for (i=start; i<=end; i++) {
237                if (get_elem_long(i)!=other.get_elem_short(i)) {
238                    *firstDifferentPos = i;
239                    result = 1;
240                    break;
241                }
242            }
243        }
244        else {
245            for (i=start; i<=end; i++) {
246                if (get_elem_long(i)!=other.get_elem_long(i)) {
247                    *firstDifferentPos = i;
248                    result = 1;
249                    break;
250                }
251            }
252        }
253    }
254
255    return result;
256}
257int ED4_bases_table::lastDifference(const ED4_bases_table& other, int start, int end, int *lastDifferentPos) const
258{
259    int i;
260    int result = 0;
261
262    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
263        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
264            for (i=end; i>=start; i--) {
265                if (get_elem_short(i)!=other.get_elem_short(i)) {
266                    *lastDifferentPos = i;
267                    result = 1;
268                    break;
269                }
270            }
271        }
272        else {
273            for (i=end; i>=start; i--) {
274                if (get_elem_short(i)!=other.get_elem_long(i)) {
275                    *lastDifferentPos = i;
276                    result = 1;
277                    break;
278                }
279            }
280        }
281    }
282    else {
283        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
284            for (i=end; i>=start; i--) {
285                if (get_elem_long(i)!=other.get_elem_short(i)) {
286                    *lastDifferentPos = i;
287                    result = 1;
288                    break;
289                }
290            }
291        }
292        else {
293            for (i=end; i>=start; i--) {
294                if (get_elem_long(i)!=other.get_elem_long(i)) {
295                    *lastDifferentPos = i;
296                    result = 1;
297                    break;
298                }
299            }
300        }
301    }
302
303    return result;
304}
305
306
307ED4_bases_table::ED4_bases_table(int maxseqlength)
308{
309    memset((char*)this, 0, sizeof(*this));
310    table_entry_size = SHORT_TABLE_ELEM_SIZE;
311    if (maxseqlength) init(maxseqlength);
312}
313
314ED4_bases_table::~ED4_bases_table()
315{
316    delete [] no_of_bases.shortTable;
317#ifdef COUNT_BASES_TABLE_SIZE
318    bases_table_size -= (no_of_entries+1)*table_entry_size;
319    printf("bases_table_size = %li\n", bases_table_size);
320#endif
321}
322
323void ED4_bases_table::change_table_length(int new_length, int default_entry)
324{
325    if (new_length!=no_of_entries) {
326        int min_length = new_length<no_of_entries ? new_length : no_of_entries;
327        int growth = new_length-no_of_entries;
328
329        if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
330            unsigned char *new_table = new unsigned char[new_length+1];
331
332            e4_assert(default_entry>=0 && default_entry<256);
333
334            memcpy(new_table, no_of_bases.shortTable, min_length*sizeof(*new_table));
335            new_table[new_length] = no_of_bases.shortTable[no_of_entries];
336            if (growth>0) {
337                for (int e=no_of_entries; e<new_length; ++e) new_table[e] = default_entry;
338            }
339
340            delete [] no_of_bases.shortTable;
341            no_of_bases.shortTable = new_table;
342        }
343        else {
344            int *new_table = new int[new_length+1];
345
346            memcpy(new_table, no_of_bases.longTable, min_length*sizeof(*new_table));
347            new_table[new_length] = no_of_bases.longTable[no_of_entries];
348            if (growth>0) {
349                for (int e=no_of_entries; e<new_length; ++e) new_table[e] = default_entry;
350            }
351
352            delete [] no_of_bases.longTable;
353            no_of_bases.longTable = new_table;
354        }
355
356#ifdef COUNT_BASES_TABLE_SIZE
357        bases_table_size += growth*table_entry_size;
358        printf("bases_table_size = %li\n", bases_table_size);
359#endif
360        no_of_entries = new_length;
361    }
362}
363
364#if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED)
365
366int ED4_bases_table::empty() const
367{
368    int i;
369
370    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
371        for (i=0; i<no_of_entries; i++) {
372            if (get_elem_short(i)) {
373                return 0;
374            }
375        }
376    }
377    else {
378        for (i=0; i<no_of_entries; i++) {
379            if (get_elem_long(i)) {
380                return 0;
381            }
382        }
383    }
384
385    return 1;
386}
387#endif // ASSERTION_USED
388
389// ------------------------
390//      Build consensus
391
392// we make static copies of the awars to avoid performance breakdown (BK_up_to_date is changed by callback ED4_consensus_definition_changed)
393
394struct ConsensusBuildParams {
395    int countgaps;
396    int gapbound;
397    int group;
398    int considbound;
399    int upper;
400    int lower;
401
402    ConsensusBuildParams(AW_root *awr) {
403        countgaps   = awr->awar(ED4_AWAR_CONSENSUS_COUNTGAPS)->read_int();
404        gapbound    = awr->awar(ED4_AWAR_CONSENSUS_GAPBOUND)->read_int();
405        group       = awr->awar(ED4_AWAR_CONSENSUS_GROUP)->read_int();
406        considbound = awr->awar(ED4_AWAR_CONSENSUS_CONSIDBOUND)->read_int();
407        upper       = awr->awar(ED4_AWAR_CONSENSUS_UPPER)->read_int();
408        lower       = awr->awar(ED4_AWAR_CONSENSUS_LOWER)->read_int();
409    }
410#if defined(UNIT_TESTS) // UT_DIFF
411    ConsensusBuildParams() {
412        // uses awar defaults
413        countgaps   = 1;
414        gapbound    = 60;
415        group       = 1;
416        considbound = 30;
417        upper       = 95;
418        lower       = 70;
419    }
420#endif
421};
422
423static ConsensusBuildParams *BK = NULL; // @@@ make member of ED4_char_table ?
424
425void ED4_consensus_definition_changed(AW_root*) {
426    delete BK; BK = 0; // invalidate
427    ED4_ROOT->request_refresh_for_consensus_terminals();
428}
429
430static ARB_ERROR toggle_consensus_display(ED4_base *base, AW_CL show) {
431    if (base->is_consensus_manager()) {
432        ED4_manager *consensus_man = base->to_manager();
433        ED4_spacer_terminal *spacer = consensus_man->parent->get_defined_level(ED4_L_SPACER)->to_spacer_terminal();
434
435        if (show) {
436            consensus_man->unhide_children();
437            spacer->extension.size[HEIGHT] = SPACERHEIGHT;
438        }
439        else {
440            consensus_man->hide_children();
441
442            ED4_group_manager *group_man = consensus_man->get_parent(ED4_L_GROUP)->to_group_manager();
443            spacer->extension.size[HEIGHT] = (group_man->dynamic_prop&ED4_P_IS_FOLDED) ? SPACERNOCONSENSUSHEIGHT : SPACERHEIGHT;
444        }
445    }
446
447    return NULL;
448}
449
450void ED4_consensus_display_changed(AW_root *root) {
451    int show = root->awar(ED4_AWAR_CONSENSUS_SHOW)->read_int();
452    ED4_ROOT->root_group_man->route_down_hierarchy(toggle_consensus_display, show).expect_no_error();
453}
454
455char *ED4_char_table::build_consensus_string(PosRange r) const {
456    ExplicitRange  range(r, size());
457    long           entries = range.size();
458    char          *new_buf = (char*)malloc(entries+1);
459
460    build_consensus_string_to(new_buf, range);
461    new_buf[entries] = 0;
462
463    return new_buf;
464}
465
466void ED4_char_table::build_consensus_string_to(char *consensus_string, ExplicitRange range) const {
467    // 'consensus_string' has to be a buffer of size 'range.size()+1'
468    // Note : Always check that consensus behavior is identical to that used in CON_evaluatestatistic()
469
470    if (!BK) BK = new ConsensusBuildParams(ED4_ROOT->aw_root);
471
472    e4_assert(consensus_string);
473    e4_assert(range.end()<size());
474
475#define PERCENT(part, all)      ((100*(part))/(all))
476#define MAX_BASES_TABLES 41     // 25
477
478    e4_assert(used_bases_tables<=MAX_BASES_TABLES);     // this is correct for DNA/RNA -> build_consensus_string() must be corrected for AMI/PRO
479
480    const int left_idx  = range.start();
481    const int right_idx = range.end();
482
483    if (sequences) {
484        for (int i=left_idx; i<=right_idx; i++) {
485            int o        = i-left_idx;
486            int bases    = 0;           // count of all bases together
487            int base[MAX_BASES_TABLES]; // count of specific base
488            int j;
489            int max_base = -1;          // maximum count of specific base
490            int max_j    = -1;          // index of this base
491
492            for (j=0; j<used_bases_tables; j++) {
493                base[j] = linear_table(j)[i];
494                if (!ADPP_IS_ALIGN_CHARACTER(index_to_upperChar(j))) {
495                    bases += base[j];
496                    if (base[j]>max_base) { // search for most used base
497                        max_base = base[j];
498                        max_j = j;
499                    }
500                }
501            }
502
503            int gaps = sequences-bases; // count of gaps
504
505            // What to do with gaps?
506
507            if (gaps == sequences) {
508                consensus_string[o] = '=';
509            }
510            else if (BK->countgaps && PERCENT(gaps, sequences)>=BK->gapbound) {
511                consensus_string[o] = '-';
512            }
513            else {
514                // Simplify using IUPAC :
515
516                char kchar;     // character for consensus
517                int  kcount;    // count this character
518
519                if (BK->group) { // group -> iupac
520                    if (ali_type == GB_AT_RNA || ali_type == GB_AT_DNA) {
521                        int no_of_bases = 0; // count of different bases used to create iupac
522                        char used_bases[MAX_BASES_TABLES+1]; // string containing those bases
523
524                        kcount = 0;
525                        for (j=0; j<used_bases_tables; j++) {
526                            int bchar = index_to_upperChar(j);
527
528                            if (!ADPP_IS_ALIGN_CHARACTER(bchar)) {
529                                if (PERCENT(base[j], sequences)>=BK->considbound) {
530                                    used_bases[no_of_bases++] = index_to_upperChar(j);
531                                    kcount += base[j];
532                                }
533                            }
534                        }
535                        used_bases[no_of_bases] = 0;
536                        kchar = ED4_encode_iupac(used_bases, ali_type);
537                    }
538                    else {
539                        e4_assert(ali_type == GB_AT_AA);
540
541                        const int amino_groups = iupac::AA_GROUP_COUNT;
542                        int       group_count[amino_groups];
543
544                        for (j=0; j<amino_groups; j++) {
545                            group_count[j] = 0;
546                        }
547                        for (j=0; j<used_bases_tables; j++) {
548                            unsigned char bchar = index_to_upperChar(j);
549
550                            if (!ADPP_IS_ALIGN_CHARACTER(bchar)) {
551                                if (PERCENT(base[j], sequences)>=BK->considbound) {
552                                    group_count[iupac::get_amino_group_for(bchar)] += base[j];
553                                }
554                            }
555                        }
556
557                        kcount = 0;
558                        int bestGroup = 0;
559                        for (j=0; j<amino_groups; j++) {
560                            if (group_count[j]>kcount) {
561                                bestGroup = j;
562                                kcount = group_count[j];
563                            }
564                        }
565
566                        kchar = iupac::get_amino_consensus_char(iupac::Amino_Group(bestGroup));
567                    }
568                }
569                else {
570                    e4_assert(max_base);            // expect at least one base to occur
571                    e4_assert(max_j >= 0);
572
573                    kchar  = index_to_upperChar(max_j);
574                    kcount = max_base;
575                    e4_assert(kchar);
576                }
577
578                // show as upper or lower case ?
579
580                int percent = PERCENT(kcount, sequences);
581
582                if (percent>=BK->upper) {
583                    consensus_string[o] = kchar;
584                }
585                else if (percent>=BK->lower) {
586                    consensus_string[o] = tolower(kchar);
587                }
588                else {
589                    consensus_string[o] = '.';
590                }
591            }
592            e4_assert(consensus_string[o]);
593        }
594    }
595    else {
596        memset(consensus_string, '?', right_idx-left_idx+1);
597    }
598#undef PERCENT
599}
600
601// -----------------------
602//      ED4_char_table
603
604bool               ED4_char_table::initialized       = false;
605unsigned char      ED4_char_table::char_to_index_tab[MAXCHARTABLE];
606unsigned char     *ED4_char_table::upper_index_chars = 0;
607unsigned char     *ED4_char_table::lower_index_chars = 0;
608int                ED4_char_table::used_bases_tables = 0;
609GB_alignment_type  ED4_char_table::ali_type          = GB_AT_RNA;
610
611inline void ED4_char_table::set_char_to_index(unsigned char c, int index)
612{
613    e4_assert(index>=0 && index<used_bases_tables);
614    char_to_index_tab[upper_index_chars[index] = toupper(c)] = index;
615    char_to_index_tab[lower_index_chars[index] = tolower(c)] = index;
616}
617unsigned char ED4_char_table::index_to_upperChar(int index) const
618{
619    return upper_index_chars[index];
620}
621unsigned char ED4_char_table::index_to_lowerChar(int index) const
622{
623    return lower_index_chars[index];
624}
625
626void ED4_char_table::expand_tables() {
627    int i;
628    for (i=0; i<used_bases_tables; i++) {
629        linear_table(i).expand_table_entry_size();
630    }
631}
632
633void ED4_char_table::initial_setup(const char *gap_chars, GB_alignment_type ali_type_) {
634    const char *groups = 0;
635    used_bases_tables  = strlen(gap_chars);
636    ali_type           = ali_type_;
637
638    switch (ali_type) {
639        case GB_AT_RNA:
640            groups = "A,C,G,TU,MRWSYKVHDBN,";
641            break;
642        case GB_AT_DNA:
643            groups = "A,C,G,UT,MRWSYKVHDBN,";
644            break;
645        case GB_AT_AA:
646            groups = "P,A,G,S,T,Q,N,E,D,B,Z,H,K,R,L,I,V,M,F,Y,W,C,X"; // @@@ DRY (create 'groups' from AA_GROUP_...)
647            break;
648        default:
649            e4_assert(0);
650            break;
651    }
652
653    e4_assert(groups);
654
655    int i;
656
657    for (i=0; groups[i]; i++) {
658        if (groups[i]==',') {
659            used_bases_tables++;
660        }
661    }
662
663    lower_index_chars = new unsigned char[used_bases_tables];
664    upper_index_chars = new unsigned char[used_bases_tables];
665
666    int idx = 0;
667    unsigned char init_val = used_bases_tables-1; // map all unknown stuff to last table
668    memset(char_to_index_tab, init_val, MAXCHARTABLE*sizeof(*char_to_index_tab));
669
670    for (i=0; groups[i]; i++) {
671        if (groups[i]==',') {
672            idx++;
673        }
674        else {
675            e4_assert(isupper(groups[i]));
676            set_char_to_index(groups[i], idx);
677        }
678    }
679
680    const char *align_string_ptr = gap_chars;
681    while (1) {
682        char c = *align_string_ptr++;
683        if (!c) break;
684        set_char_to_index(c, idx++);
685    }
686
687    e4_assert(idx==used_bases_tables);
688    initialized = true;
689}
690
691ED4_char_table::ED4_char_table(int maxseqlength)
692{
693    memset((char*)this, 0, sizeof(*this));
694
695    if (!initialized) {
696        char *align_string = ED4_ROOT->aw_root->awar_string(ED4_AWAR_GAP_CHARS)->read_string();
697        initial_setup(align_string, ED4_ROOT->alignment_type);
698        free(align_string);
699    }
700
701    bases_table = new ED4_bases_table_ptr[used_bases_tables];
702
703    int i;
704    for (i=0; i<used_bases_tables; i++) {
705        bases_table[i] = new ED4_bases_table(maxseqlength);
706    }
707
708    sequences = 0;
709}
710
711void ED4_char_table::init(int maxseqlength)
712{
713    int i;
714    for (i=0; i<used_bases_tables; i++) {
715        linear_table(i).init(maxseqlength);
716    }
717
718    sequences = 0;
719}
720
721void ED4_char_table::bases_and_gaps_at(int column, int *bases, int *gaps) const
722{
723    int i,
724        b = 0,
725        g = 0;
726
727    for (i=0; i<used_bases_tables; i++) {
728        char c = upper_index_chars[i];
729
730        if (ADPP_IS_ALIGN_CHARACTER(c)) {
731            g += table(c)[column];
732        }
733        else {
734            b += table(c)[column];
735        }
736    }
737
738    if (bases) *bases = b;
739    if (gaps)  *gaps  = g;
740}
741
742ED4_char_table::~ED4_char_table()
743{
744    int i;
745    for (i=0; i<used_bases_tables; i++) {
746        delete bases_table[i];
747    }
748
749    delete [] bases_table;
750}
751
752const PosRange *ED4_char_table::changed_range(const ED4_char_table& other) const
753{
754    int i;
755    int Size = size();
756    int start = Size-1;
757    int end = 0;
758
759    e4_assert(Size==other.size());
760    for (i=0; i<used_bases_tables; i++) {
761        if (linear_table(i).firstDifference(other.linear_table(i), 0, start, &start)) {
762            if (end<start) {
763                end = start;
764            }
765            linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end);
766
767            for (; i<used_bases_tables; i++) {
768                linear_table(i).firstDifference(other.linear_table(i), 0, start-1, &start);
769                if (end<start) {
770                    end = start;
771                }
772                linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end);
773            }
774
775            e4_assert(start<=end);
776
777            static PosRange range;
778            range = PosRange(start, end);
779
780            return &range;
781        }
782    }
783    return NULL;
784}
785void ED4_char_table::add(const ED4_char_table& other)
786{
787    if (other.ignore) return;
788    if (other.size()==0) return;
789    prepare_to_add_elements(other.sequences);
790    add(other, 0, other.size()-1);
791}
792void ED4_char_table::add(const ED4_char_table& other, int start, int end)
793{
794    if (other.ignore) return;
795
796    test();
797    other.test();
798
799    e4_assert(start<=end);
800
801    int i;
802    for (i=0; i<used_bases_tables; i++) {
803        linear_table(i).add(other.linear_table(i), start, end);
804    }
805
806    sequences += other.sequences;
807
808    test();
809}
810
811void ED4_char_table::sub(const ED4_char_table& other)
812{
813    if (other.ignore) return;
814    sub(other, 0, other.size()-1);
815}
816
817void ED4_char_table::sub(const ED4_char_table& other, int start, int end)
818{
819    if (other.ignore) return;
820
821    test();
822    other.test();
823    e4_assert(start<=end);
824
825    int i;
826    for (i=0; i<used_bases_tables; i++) {
827        linear_table(i).sub(other.linear_table(i), start, end);
828    }
829
830    sequences -= other.sequences;
831
832    test();
833}
834
835void ED4_char_table::sub_and_add(const ED4_char_table& Sub, const ED4_char_table& Add) {
836    test();
837
838    if (Sub.ignore) {
839        e4_assert(Add.ignore);
840        return;
841    }
842
843    Sub.test();
844    Add.test();
845
846    const PosRange *range = Sub.changed_range(Add);
847    if (range) {
848        prepare_to_add_elements(Add.added_sequences()-Sub.added_sequences());
849        sub_and_add(Sub, Add, *range);
850    }
851
852    test();
853}
854
855void ED4_char_table::sub_and_add(const ED4_char_table& Sub, const ED4_char_table& Add, PosRange range) {
856    test();
857    Sub.test();
858    Add.test();
859
860    e4_assert(!Sub.ignore && !Add.ignore);
861    e4_assert(range.is_part());
862
863    int i;
864    for (i=0; i<used_bases_tables; i++) {
865        linear_table(i).sub_and_add(Sub.linear_table(i), Add.linear_table(i), range);
866    }
867
868    test();
869}
870
871const PosRange *ED4_char_table::changed_range(const char *string1, const char *string2, int min_len)
872{
873    const unsigned long *range1 = (const unsigned long *)string1;
874    const unsigned long *range2 = (const unsigned long *)string2;
875    const int            step   = sizeof(*range1);
876
877    int l = min_len/step;
878    int r = min_len%step;
879    int i;
880    int j;
881    int k;
882
883    int start = -1, end = -1;
884
885    for (i=0; i<l; i++) {       // search wordwise (it's faster)
886        if (range1[i] != range2[i]) {
887            k = i*step;
888            for (j=0; j<step; j++) {
889                if (string1[k+j] != string2[k+j]) {
890                    start = end = k+j;
891                    break;
892                }
893            }
894
895            break;
896        }
897    }
898
899    k = l*step;
900    if (i==l) {         // no difference found in word -> look at rest
901        for (j=0; j<r; j++) {
902            if (string1[k+j] != string2[k+j]) {
903                start = end = k+j;
904                break;
905            }
906        }
907
908        if (j==r) {     // the strings are equal
909            return 0;
910        }
911    }
912
913    e4_assert(start != -1 && end != -1);
914
915    for (j=r-1; j>=0; j--) {    // search rest backwards
916        if (string1[k+j] != string2[k+j]) {
917            end = k+j;
918            break;
919        }
920    }
921
922    if (j==-1) {                // not found in rest -> search words backwards
923        int m = start/step;
924        for (i=l-1; i>=m; i--) {
925            if (range1[i] != range2[i]) {
926                k = i*step;
927                for (j=step-1; j>=0; j--) {
928                    if (string1[k+j] != string2[k+j]) {
929                        end = k+j;
930                        break;
931                    }
932                }
933                break;
934            }
935        }
936    }
937
938    e4_assert(start<=end);
939
940    static PosRange range;
941    range = PosRange(start, end);
942
943    return &range;
944}
945
946void ED4_char_table::add(const char *scan_string, int len)
947{
948    test();
949
950    prepare_to_add_elements(1);
951
952    int i;
953    int sz = size();
954
955    if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) {
956        for (i=0; i<len; i++) {
957            unsigned char c = scan_string[i];
958            if (!c) break;
959            table(c).inc_short(i);
960        }
961        ED4_bases_table& t = table('.');
962        for (; i<sz; i++) {
963            t.inc_short(i);
964        }
965    }
966    else {
967        for (i=0; i<len; i++) {
968            unsigned char c = scan_string[i];
969            if (!c) break;
970            table(c).inc_long(i);
971        }
972        ED4_bases_table& t = table('.');
973        for (; i<sz; i++) {
974            t.inc_long(i);
975        }
976    }
977
978    sequences++;
979
980    test();
981}
982
983void ED4_char_table::sub(const char *scan_string, int len)
984{
985    test();
986
987    int i;
988    int sz = size();
989
990    if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) {
991        for (i=0; i<len; i++) {
992            unsigned char c = scan_string[i];
993            if (!c) break;
994            table(c).dec_short(i);
995        }
996        ED4_bases_table& t = table('.');
997        for (; i<sz; i++) {
998            t.dec_short(i);
999        }
1000    }
1001    else {
1002        for (i=0; i<len; i++) {
1003            unsigned char c = scan_string[i];
1004            if (!c) break;
1005            table(c).dec_long(i);
1006        }
1007        ED4_bases_table& t = table('.');
1008        for (; i<sz; i++) {
1009            t.dec_long(i);
1010        }
1011    }
1012
1013    sequences--;
1014
1015    test();
1016}
1017
1018void ED4_char_table::sub_and_add(const char *old_string, const char *new_string, PosRange range) {
1019    test();
1020    int start = range.start();
1021    int end   = range.end();
1022    e4_assert(start<=end);
1023
1024    int i;
1025    ED4_bases_table& t = table('.');
1026
1027    if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) {
1028        for (i=start; i<=end; i++) {
1029            unsigned char o = old_string[i];
1030            unsigned char n = new_string[i];
1031
1032            e4_assert(o); // @@@ if these never fail, some code below is superfluous
1033            e4_assert(n);
1034           
1035            if (!o) {
1036                for (; n && i<=end; i++) {
1037                    n = new_string[i];
1038                    table(n).inc_short(i);
1039                    t.dec_short(i);
1040                }
1041            }
1042            else if (!n) {
1043                for (; o && i<=end; i++) {
1044                    o = old_string[i];
1045                    table(o).dec_short(i);
1046                    t.inc_short(i);
1047                }
1048
1049            }
1050            else if (o!=n) {
1051                table(o).dec_short(i);
1052                table(n).inc_short(i);
1053
1054            }
1055        }
1056    }
1057    else {
1058        for (i=start; i<=end; i++) {
1059            unsigned char o = old_string[i];
1060            unsigned char n = new_string[i];
1061
1062            e4_assert(o); // @@@ if these never fail, some code below is superfluous
1063            e4_assert(n);
1064           
1065            if (!o) {
1066                for (; n && i<=end; i++) {
1067                    n = new_string[i];
1068                    table(n).inc_long(i);
1069                    t.dec_long(i);
1070                }
1071            }
1072            else if (!n) {
1073                for (; o && i<=end; i++) {
1074                    o = old_string[i];
1075                    table(o).dec_long(i);
1076                    t.inc_long(i);
1077                }
1078
1079            }
1080            else if (o!=n) {
1081                table(o).dec_long(i);
1082                table(n).inc_long(i);
1083            }
1084        }
1085    }
1086
1087    test();
1088}
1089
1090void ED4_char_table::change_table_length(int new_length) {
1091    for (int c=0; c<used_bases_tables; c++) {
1092        linear_table(c).change_table_length(new_length, 0);
1093    }
1094    test();
1095}
1096
1097//  --------------
1098//      tests:
1099
1100#if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED)
1101
1102bool ED4_char_table::empty() const
1103{
1104    int i;
1105    for (i=0; i<used_bases_tables; i++) {
1106        if (!linear_table(i).empty()) {
1107            return false;
1108        }
1109    }
1110    return true;
1111}
1112
1113bool ED4_char_table::ok() const
1114{
1115    if (empty()) return true;
1116    if (is_ignored()) return true;
1117
1118    if (sequences<0) {
1119        fprintf(stderr, "Negative # of sequences (%i) in ED4_char_table\n", sequences);
1120        return false;
1121    }
1122
1123    int i;
1124    for (i=0; i<MAXSEQUENCECHARACTERLENGTH; i++) {
1125        int bases;
1126        int gaps;
1127
1128        bases_and_gaps_at(i, &bases, &gaps);
1129
1130        if (!bases && !gaps) { // this occurs after insert column
1131            int j;
1132            for (j=i+1; j<MAXSEQUENCECHARACTERLENGTH; j++) {    // test if we find any bases till end of table !!!
1133
1134                bases_and_gaps_at(j, &bases, 0);
1135                if (bases) { // bases found -> empty position was illegal
1136                    fprintf(stderr, "Zero in ED4_char_table at column %i\n", i);
1137                    return false;
1138                }
1139            }
1140            break;
1141        }
1142        else {
1143            if ((bases+gaps)!=sequences) {
1144                fprintf(stderr, "bases+gaps should be equal to # of sequences at column %i\n", i);
1145                return false;
1146            }
1147        }
1148    }
1149
1150    return true;
1151}
1152
1153#endif // TEST_CHAR_TABLE_INTEGRITY || ASSERTION_USED
1154
1155#if defined(TEST_CHAR_TABLE_INTEGRITY)
1156
1157void ED4_char_table::test() const {
1158
1159    if (!ok()) {
1160        GBK_terminate("ED4_char_table::test() failed");
1161    }
1162}
1163
1164#endif // TEST_CHAR_TABLE_INTEGRITY
1165
1166
1167// --------------------------------------------------------------------------------
1168
1169#ifdef UNIT_TESTS
1170#ifndef TEST_UNIT_H
1171#include <test_unit.h>
1172#endif
1173
1174void TEST_char_table() {
1175    const char alphabeth[]   = "ACGTN-";
1176    const int alphabeth_size = strlen(alphabeth);
1177
1178    const int seqlen = 30;
1179    char      seq[seqlen+1];
1180
1181    const char *gapChars = ".-"; // @@@ doesnt make any difference which gaps are used - why?
1182    // const char *gapChars = "-";
1183    // const char *gapChars = ".";
1184    // const char *gapChars = "A"; // makes me fail
1185    ED4_char_table::initial_setup(gapChars, GB_AT_RNA);
1186    ED4_init_is_align_character(gapChars);
1187
1188    TEST_REJECT(BK);
1189    BK = new ConsensusBuildParams;
1190
1191    // BK->lower = 70; BK->upper = 95; BK->gapbound = 60; // defaults from awars
1192    BK->lower    = 40; BK->upper = 70; BK->gapbound = 40;
1193
1194    srand(100);
1195    for (int loop = 0; loop<5; ++loop) {
1196        unsigned seed = rand();
1197
1198        srand(seed);
1199
1200        ED4_char_table tab(seqlen);
1201
1202        // build some seq
1203        for (int c = 0; c<seqlen; ++c) seq[c] = alphabeth[rand()%alphabeth_size];
1204        seq[seqlen]                           = 0;
1205
1206        TEST_EXPECT_EQUAL(strlen(seq), size_t(seqlen));
1207
1208        const int  add_count = 300;
1209        char      *added_seqs[add_count];
1210
1211        // add seq multiple times
1212        for (int a = 0; a<add_count; ++a) {
1213            tab.add(seq, seqlen);
1214            added_seqs[a] = strdup(seq);
1215
1216            // modify 1 character in seq:
1217            int sidx  = rand()%seqlen;
1218            int aidx  = rand()%alphabeth_size;
1219            seq[sidx] = alphabeth[aidx];
1220        }
1221
1222        // build consensi (just check regression)
1223        {
1224            char *consensus = tab.build_consensus_string();
1225            switch (seed) {
1226                case 677741240: TEST_EXPECT_EQUAL(consensus, ".-s-NW..aWu.NnWYa.R.mgcNK.c..."); break;
1227                case 721151648: TEST_EXPECT_EQUAL(consensus, "a.nn..K..-gU.RW-SNcau.WNNacn.u"); break;
1228                case 345295160: TEST_EXPECT_EQUAL(consensus, "..-g...MSn...guc.n.u.R.n.-Ng.k"); break;
1229                case 346389111: TEST_EXPECT_EQUAL(consensus, ".unAn...gN.kc-cS.Raun...Sa-gY."); break;
1230                case 367171911: TEST_EXPECT_EQUAL(consensus, "na.Na.nu.c-.-NU.aYgn-nng-.Wa.M"); break;
1231                default: TEST_EXPECT_EQUAL(consensus, "undef");
1232            }
1233            free(consensus);
1234        }
1235
1236        for (int a = 0; a<add_count; a++) {
1237            tab.sub(added_seqs[a], seqlen);
1238            freenull(added_seqs[a]);
1239        }
1240        {
1241            char *consensus = tab.build_consensus_string();
1242            TEST_EXPECT_EQUAL(consensus, "??????????????????????????????"); // check tab is empty
1243            free(consensus);
1244        }
1245    }
1246
1247    delete BK;
1248    BK = 0;
1249}
1250
1251#endif // UNIT_TESTS
1252
1253// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.