source: branches/alilink/SL/CONSENSUS/chartable.cxx

Last change on this file was 18126, checked in by westram, 5 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 44.3 KB
Line 
1// ================================================================= //
2//                                                                   //
3//   File      : chartable.cxx                                       //
4//   Purpose   : count characters of multiple sequences              //
5//                                                                   //
6//   Institute of Microbiology (Technical University Munich)         //
7//   http://www.arb-home.de/                                         //
8//                                                                   //
9// ================================================================= //
10
11#include "chartable.h"
12
13#include <arb_msg.h>
14#include <arb_mem.h>
15#include <iupac.h>
16#include <arb_defs.h>
17#include <smartptr.h>
18
19#include <cctype>
20
21#define CONSENSUS_AWAR_SOURCE CAS_INTERNAL // not awar-constructable here
22#include <consensus.h>
23
24using namespace chartable;
25
26// ------------------------
27//      SepBaseFreq
28
29#ifdef DEBUG
30// # define COUNT_BASES_TABLE_SIZE
31#endif
32
33#ifdef COUNT_BASES_TABLE_SIZE
34static long bases_table_size = 0L;
35#endif
36
37void SepBaseFreq::init(int length) {
38    if (length) {
39        if (no_of_bases.shortTable) {
40            delete no_of_bases.shortTable;
41            no_of_bases.shortTable = NULp;
42#ifdef COUNT_BASES_TABLE_SIZE
43            bases_table_size -= (no_of_entries+1)*table_entry_size;
44#endif
45        }
46
47        no_of_entries = length;
48        int new_size = (no_of_entries+1)*table_entry_size;
49        no_of_bases.shortTable = (unsigned char*)new char[new_size];
50        memset(no_of_bases.shortTable, 0, new_size);
51#ifdef COUNT_BASES_TABLE_SIZE
52        bases_table_size += new_size;
53        printf("bases_table_size = %li\n", bases_table_size);
54#endif
55    }
56}
57
58void SepBaseFreq::expand_table_entry_size() { // converts short table into long table
59    ct_assert(table_entry_size==SHORT_TABLE_ELEM_SIZE);
60
61    int *new_table = new int[no_of_entries+1];
62    int i;
63
64    for (i=0; i<=no_of_entries; i++) { // LOOP_VECTORIZED[!<5.0]
65        new_table[i] = no_of_bases.shortTable[i];
66    }
67    delete [] no_of_bases.shortTable;
68
69#ifdef COUNT_BASES_TABLE_SIZE
70    bases_table_size += (no_of_entries+1)*(LONG_TABLE_ELEM_SIZE-SHORT_TABLE_ELEM_SIZE);
71    printf("bases_table_size = %li\n", bases_table_size);
72#endif
73
74    no_of_bases.longTable = new_table;
75    table_entry_size = LONG_TABLE_ELEM_SIZE;
76}
77
78#define INVALID_TABLE_OPERATION() GBK_terminatef("SepBaseFreq: invalid operation at %i", __LINE__)
79
80void SepBaseFreq::add(const SepBaseFreq& other, int start, int end) {
81    ct_assert(no_of_entries==other.no_of_entries);
82    ct_assert(start>=0);
83    ct_assert(end<no_of_entries);
84    ct_assert(start<=end);
85
86    int i;
87    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
88        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
89            for (i=start; i<=end; i++) {
90                set_elem_short(i, get_elem_short(i)+other.get_elem_short(i));
91            }
92        }
93        else {
94            INVALID_TABLE_OPERATION(); // cannot add long to short
95        }
96    }
97    else {
98        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
99            for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
100                set_elem_long(i, get_elem_long(i)+other.get_elem_short(i));
101            }
102        }
103        else {
104            for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
105                set_elem_long(i, get_elem_long(i)+other.get_elem_long(i));
106            }
107        }
108    }
109}
110void SepBaseFreq::sub(const SepBaseFreq& other, int start, int end) {
111    ct_assert(no_of_entries==other.no_of_entries);
112    ct_assert(start>=0);
113    ct_assert(end<no_of_entries);
114    ct_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            INVALID_TABLE_OPERATION(); // cannot sub long from short
125        }
126    }
127    else {
128        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
129            for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
130                set_elem_long(i, get_elem_long(i)-other.get_elem_short(i));
131            }
132        }
133        else {
134            for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
135                set_elem_long(i, get_elem_long(i)-other.get_elem_long(i));
136            }
137        }
138    }
139}
140void SepBaseFreq::sub_and_add(const SepBaseFreq& Sub, const SepBaseFreq& Add, PosRange range) {
141    ct_assert(no_of_entries==Sub.no_of_entries);
142    ct_assert(no_of_entries==Add.no_of_entries);
143
144    int start = range.start();
145    int end   = range.end();
146
147    ct_assert(start>=0);
148    ct_assert(end<no_of_entries);
149    ct_assert(start<=end);
150
151    int i;
152    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
153        if (Sub.table_entry_size==SHORT_TABLE_ELEM_SIZE &&
154            Add.table_entry_size==SHORT_TABLE_ELEM_SIZE)
155        {
156            for (i=start; i<=end; i++) {
157                set_elem_short(i, get_elem_short(i)-Sub.get_elem_short(i)+Add.get_elem_short(i));
158            }
159        }
160        else {
161            INVALID_TABLE_OPERATION(); // cannot add or sub long to/from short
162        }
163    }
164    else {
165        if (Sub.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
166            if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
167                for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
168                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_short(i));
169                }
170            }
171            else {
172                for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
173                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_short(i)+Add.get_elem_long(i));
174                }
175            }
176        }
177        else {
178            if (Add.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
179                for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
180                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_short(i));
181                }
182            }
183            else {
184                for (i=start; i<=end; i++) { // LOOP_VECTORIZED[!>=620<7]
185                    set_elem_long(i, get_elem_long(i)-Sub.get_elem_long(i)+Add.get_elem_long(i));
186                }
187            }
188        }
189    }
190}
191
192int SepBaseFreq::firstDifference(const SepBaseFreq& other, int start, int end, int *firstDifferentPos) const {
193    int i;
194    int result = 0;
195
196    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
197        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
198            for (i=start; i<=end; i++) {
199                if (get_elem_short(i)!=other.get_elem_short(i)) {
200                    *firstDifferentPos = i;
201                    result = 1;
202                    break;
203                }
204            }
205        }
206        else {
207            for (i=start; i<=end; i++) {
208                if (get_elem_short(i)!=other.get_elem_long(i)) {
209                    *firstDifferentPos = i;
210                    result = 1;
211                    break;
212                }
213            }
214        }
215    }
216    else {
217        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
218            for (i=start; i<=end; i++) {
219                if (get_elem_long(i)!=other.get_elem_short(i)) {
220                    *firstDifferentPos = i;
221                    result = 1;
222                    break;
223                }
224            }
225        }
226        else {
227            for (i=start; i<=end; i++) {
228                if (get_elem_long(i)!=other.get_elem_long(i)) {
229                    *firstDifferentPos = i;
230                    result = 1;
231                    break;
232                }
233            }
234        }
235    }
236
237    return result;
238}
239int SepBaseFreq::lastDifference(const SepBaseFreq& other, int start, int end, int *lastDifferentPos) const {
240    int i;
241    int result = 0;
242
243    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
244        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
245            for (i=end; i>=start; i--) {
246                if (get_elem_short(i)!=other.get_elem_short(i)) {
247                    *lastDifferentPos = i;
248                    result = 1;
249                    break;
250                }
251            }
252        }
253        else {
254            for (i=end; i>=start; i--) {
255                if (get_elem_short(i)!=other.get_elem_long(i)) {
256                    *lastDifferentPos = i;
257                    result = 1;
258                    break;
259                }
260            }
261        }
262    }
263    else {
264        if (other.table_entry_size==SHORT_TABLE_ELEM_SIZE) {
265            for (i=end; i>=start; i--) {
266                if (get_elem_long(i)!=other.get_elem_short(i)) {
267                    *lastDifferentPos = i;
268                    result = 1;
269                    break;
270                }
271            }
272        }
273        else {
274            for (i=end; i>=start; i--) {
275                if (get_elem_long(i)!=other.get_elem_long(i)) {
276                    *lastDifferentPos = i;
277                    result = 1;
278                    break;
279                }
280            }
281        }
282    }
283
284    return result;
285}
286
287
288SepBaseFreq::SepBaseFreq(int maxseqlength)
289    : table_entry_size(SHORT_TABLE_ELEM_SIZE),
290      no_of_entries(0)
291{
292    no_of_bases.shortTable = NULp;
293    if (maxseqlength) init(maxseqlength);
294}
295
296SepBaseFreq::~SepBaseFreq() {
297    delete [] no_of_bases.shortTable;
298#ifdef COUNT_BASES_TABLE_SIZE
299    bases_table_size -= (no_of_entries+1)*table_entry_size;
300    printf("bases_table_size = %li\n", bases_table_size);
301#endif
302}
303
304void SepBaseFreq::change_table_length(int new_length, int default_entry) {
305    if (new_length!=no_of_entries) {
306        int min_length = new_length<no_of_entries ? new_length : no_of_entries;
307        int growth = new_length-no_of_entries;
308
309        if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
310            unsigned char *new_table = new unsigned char[new_length+1];
311
312            ct_assert(default_entry>=0 && default_entry<256);
313
314            memcpy(new_table, no_of_bases.shortTable, min_length*sizeof(*new_table));
315            new_table[new_length] = no_of_bases.shortTable[no_of_entries];
316            if (growth>0) {
317                for (int e=no_of_entries; e<new_length; ++e) new_table[e] = default_entry;
318            }
319
320            delete [] no_of_bases.shortTable;
321            no_of_bases.shortTable = new_table;
322        }
323        else {
324            int *new_table = new int[new_length+1];
325
326            memcpy(new_table, no_of_bases.longTable, min_length*sizeof(*new_table));
327            new_table[new_length] = no_of_bases.longTable[no_of_entries];
328            if (growth>0) {
329                for (int e=no_of_entries; e<new_length; ++e) { // LOOP_VECTORIZED
330                    new_table[e] = default_entry;
331                }
332            }
333
334            delete [] no_of_bases.longTable;
335            no_of_bases.longTable = new_table;
336        }
337
338#ifdef COUNT_BASES_TABLE_SIZE
339        bases_table_size += growth*table_entry_size;
340        printf("bases_table_size = %li\n", bases_table_size);
341#endif
342        no_of_entries = new_length;
343    }
344}
345
346#if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED)
347
348int SepBaseFreq::empty() const {
349    int i;
350
351    if (table_entry_size==SHORT_TABLE_ELEM_SIZE) {
352        for (i=0; i<no_of_entries; i++) {
353            if (get_elem_short(i)) {
354                return 0;
355            }
356        }
357    }
358    else {
359        for (i=0; i<no_of_entries; i++) {
360            if (get_elem_long(i)) {
361                return 0;
362            }
363        }
364    }
365
366    return 1;
367}
368#endif // ASSERTION_USED
369
370char *BaseFrequencies::build_consensus_string(PosRange r, const ConsensusBuildParams& cbp) const {
371    ExplicitRange range(r, size());
372
373    long  entries = range.size();
374    char *new_buf = ARB_alloc<char>(entries+1);
375
376    build_consensus_string_to(new_buf, range, cbp);
377    new_buf[entries] = 0;
378
379    return new_buf;
380}
381
382#if defined(DEBUG)
383// #define DEBUG_CONSENSUS
384#endif
385
386void BaseFrequencies::build_consensus_string_to(char *consensus_string, ExplicitRange range, const ConsensusBuildParams& BK) const {
387    // 'consensus_string' has to be a buffer of size 'range.size()+1'
388    // Note : Always check that consensus behavior is identical to that used in CON_evaluatestatistic()
389
390    ct_assert(consensus_string);
391    ct_assert(range.end()<size());
392
393#define PERCENT(part, all)      ((100*(part))/(all))
394#define MAX_BASES_TABLES 41     // 25
395
396    ct_assert(used_bases_tables<=MAX_BASES_TABLES);     // this is correct for DNA/RNA -> build_consensus_string() must be corrected for AMI/PRO
397
398    const int left_idx  = range.start();
399    const int right_idx = range.end();
400
401#if defined(DEBUG_CONSENSUS)
402 #define DUMPINT(var) do { if (dumpcol) fprintf(stderr, "%s=%i ", #var, var); } while(0)
403#else
404 #define DUMPINT(var)
405#endif
406
407    if (sequenceUnits) {
408        for (int i=left_idx; i<=right_idx; i++) {
409#if defined(DEBUG_CONSENSUS)
410            int  column  = i+1;
411            bool dumpcol = (column == 21);
412            DUMPINT(column);
413#endif
414
415            const int o = i-left_idx; // output offset
416
417            int bases        = 0;       // count of all bases together
418            int base[MAX_BASES_TABLES]; // count of specific base
419            int max_base     = -1;      // maximum count of specific base
420            int max_base_idx = -1;      // index of this base
421
422            for (int j=0; j<used_bases_tables; j++) {
423                base[j] = linear_table(j)[i];
424                if (!isGap(index_to_upperChar(j))) {
425                    bases += base[j];
426                    if (base[j]>max_base) { // search for most used base
427                        max_base     = base[j];
428                        max_base_idx = j;
429                    }
430                }
431            }
432
433            int gaps = sequenceUnits-bases; // count of gaps
434
435            DUMPINT(sequenceUnits);
436            DUMPINT(gaps);
437            DUMPINT(bases);
438
439            // What to do with gaps?
440
441            if (gaps == sequenceUnits) {
442                consensus_string[o] = '=';
443            }
444            else if (BK.countgaps && PERCENT(gaps, sequenceUnits) >= BK.gapbound) {
445                DUMPINT(PERCENT(gaps,sequenceUnits));
446                DUMPINT(BK.gapbound);
447                consensus_string[o] = '-';
448            }
449            else {
450                char kchar  = 0; // character for consensus
451                int  kcount = 0; // count this character
452
453                if (BK.group) { // simplify using IUPAC
454                    if (ali_type == GB_AT_RNA || ali_type == GB_AT_DNA) {
455                        int no_of_bases = 0; // count of different bases used to create iupac
456                        char used_bases[MAX_BASES_TABLES+1]; // string containing those bases
457
458                        for (int j=0; j<used_bases_tables; j++) {
459                            int bchar = index_to_upperChar(j);
460
461                            if (!isGap(bchar)) {
462                                if (base[j]>0 && PERCENT(base[j],bases) >= BK.considbound) {
463#if defined(DEBUG_CONSENSUS)
464                                    if (!kcount) DUMPINT(BK.considbound);
465#endif
466                                    used_bases[no_of_bases++]  = index_to_upperChar(j);
467                                    kcount                    += base[j];
468
469                                    DUMPINT(base[j]);
470                                    DUMPINT(PERCENT(base[j],bases));
471                                    DUMPINT(kcount);
472                                }
473                            }
474                        }
475                        used_bases[no_of_bases] = 0;
476                        kchar = iupac::encode(used_bases, ali_type);
477                    }
478                    else {
479                        ct_assert(ali_type == GB_AT_AA);
480
481                        const int amino_groups = iupac::AA_GROUP_COUNT;
482                        int       group_count[amino_groups];
483
484                        for (int j=0; j<amino_groups; j++) {
485                            group_count[j] = 0;
486                        }
487                        for (int j=0; j<used_bases_tables; j++) {
488                            unsigned char bchar = index_to_upperChar(j);
489
490                            if (!isGap(bchar)) {
491                                group_count[iupac::get_amino_group_for(bchar)] += base[j];
492                            }
493                        }
494
495                        int bestGroup = 0;
496                        for (int j=0; j<amino_groups; j++) {
497                            if (group_count[j]>kcount) {
498                                bestGroup = j;
499                                kcount = group_count[j];
500                            }
501                        }
502
503                        ct_assert(kcount>0);
504                        if (PERCENT(kcount, bases) >= BK.considbound) {
505                            kchar = iupac::get_amino_consensus_char(iupac::Amino_Group(bestGroup));
506                        }
507                        else {
508                            kchar  = 'X';
509                            kcount = bases;
510                        }
511                    }
512                }
513                if (!kcount) {           // IUPAC grouping is either off OR didnt consider any bases
514                    ct_assert(max_base); // expect at least one base to occur
515                    ct_assert(max_base_idx >= 0);
516
517                    kchar  = index_to_upperChar(max_base_idx);
518                    kcount = max_base;
519                }
520
521                ct_assert(kchar);
522                ct_assert(kcount);
523                ct_assert(kcount<=bases);
524
525                // show as upper or lower case ?
526                int percent = PERCENT(kcount, BK.countgaps ? sequenceUnits : bases);
527                DUMPINT(percent);
528                DUMPINT(BK.upper);
529                DUMPINT(BK.lower);
530                if (percent>=BK.upper) {
531                    consensus_string[o] = kchar;
532                }
533                else if (percent>=BK.lower) {
534                    consensus_string[o] = tolower(kchar);
535                }
536                else {
537                    consensus_string[o] = '.';
538                }
539            }
540            ct_assert(consensus_string[o]);
541
542#if defined(DEBUG_CONSENSUS)
543            if (dumpcol) fprintf(stderr, "result='%c'\n", consensus_string[o]);
544#endif
545        }
546    }
547    else {
548        memset(consensus_string, '?', right_idx-left_idx+1);
549    }
550#undef PERCENT
551}
552
553// -----------------------
554//      BaseFrequencies
555
556bool              BaseFrequencies::initialized       = false;
557Ambiguity         BaseFrequencies::ambiguity_table[MAX_AMBIGUITY_CODES];
558uint8_t           BaseFrequencies::unitsPerSequence  = 1;
559unsigned char     BaseFrequencies::char_to_index_tab[MAXCHARTABLE];
560bool              BaseFrequencies::is_gap[MAXCHARTABLE];
561unsigned char     BaseFrequencies::upper_index_chars[MAX_USED_BASES_TABLES];
562int               BaseFrequencies::used_bases_tables = 0;
563GB_alignment_type BaseFrequencies::ali_type          = GB_AT_RNA;
564
565inline void BaseFrequencies::set_char_to_index(unsigned char c, int index) {
566    ct_assert(index>=0 && index<used_bases_tables);
567    char_to_index_tab[upper_index_chars[index] = toupper(c)] = index;
568    char_to_index_tab[tolower(c)] = index;
569}
570
571void BaseFrequencies::expand_tables() {
572    int i;
573    for (i=0; i<used_bases_tables; i++) {
574        linear_table(i).expand_table_entry_size();
575    }
576}
577
578static const char *getAAgroupDefinition() {
579    using namespace iupac;
580    static SmartCharPtr aadef;
581
582    if (aadef.isNull()) {
583        char buffer[2*26+1];
584        int  off = 0;
585        for (int c = 'A'; c<='Z'; ++c) {
586            switch (get_amino_group_for(c)) {
587                case AA_GROUP_NONE: break; // cannot add: would add invalid characters -> add 'X' manually below.
588                case AA_GROUP_ILLEGAL: ct_assert(0); break;
589                default:
590                    buffer[off++] = c;
591                    buffer[off++] = ',';
592                    break;
593            }
594        }
595        buffer[off] = 'X';
596        buffer[off+1] = 0;
597
598        aadef = strdup(buffer);
599    }
600
601    return &*aadef;
602}
603
604void BaseFrequencies::setup(const char *gap_chars, GB_alignment_type ali_type_) {
605    const char *groups    = NULp;
606    const char *ambiguous = NULp;
607
608    ali_type = ali_type_;
609
610    switch (ali_type) {
611        case GB_AT_RNA:
612            groups           = "A,C,G,TU,"; // Note: terminal ',' defines a group for unknown characters
613            ambiguous        = "MRWSYKVHDBN";
614            unitsPerSequence = 12;
615            break;
616
617        case GB_AT_DNA:
618            groups           = "A,C,G,UT,";
619            ambiguous        = "MRWSYKVHDBN";
620            unitsPerSequence = 12;
621            break;
622
623        case GB_AT_AA:
624            groups           = getAAgroupDefinition();
625            unitsPerSequence = 1;
626            break;
627
628        default:
629            ct_assert(0);
630            break;
631    }
632
633    ct_assert(groups);
634#if defined(ASSERTION_USED)
635    if (ali_type == GB_AT_AA) {
636        ct_assert(!ambiguous);
637        // ct_assert(strlen(ambiguous) <= MAX_AMBIGUITY_CODES); // @@@ later
638    }
639    else {
640        ct_assert(strlen(ambiguous) == MAX_AMBIGUITY_CODES);
641    }
642#endif
643
644    for (int i = 0; i<MAXCHARTABLE; ++i) is_gap[i] = false;
645
646    int   gapcharcount     = 0;
647    char *unique_gap_chars = strdup(gap_chars);
648    for (int i = 0; gap_chars[i]; ++i) {
649        if (!is_gap[safeCharIndex(gap_chars[i])]) { // ignore duplicate occurrences in 'gap_chars'
650            is_gap[safeCharIndex(gap_chars[i])] = true;
651            unique_gap_chars[gapcharcount++]    = gap_chars[i];
652        }
653    }
654
655    used_bases_tables = gapcharcount;
656
657    for (int i=0; groups[i]; i++) {
658        if (groups[i]==',') {
659            used_bases_tables++;
660        }
661    }
662
663    ct_assert(used_bases_tables<=MAX_USED_BASES_TABLES);
664
665    int idx = 0;
666    unsigned char init_val = used_bases_tables-1; // map all unknown stuff to last table
667    memset(char_to_index_tab, init_val, MAXCHARTABLE*sizeof(*char_to_index_tab));
668
669    for (int i=0; groups[i]; i++) {
670        if (groups[i]==',') {
671            idx++;
672        }
673        else {
674            ct_assert(isupper(groups[i]));
675            set_char_to_index(groups[i], idx);
676        }
677    }
678
679    const char *align_string_ptr = unique_gap_chars;
680    while (1) {
681        char c = *align_string_ptr++;
682        if (!c) break;
683        set_char_to_index(c, idx++);
684    }
685
686    if (ambiguous) {
687        ct_assert(ali_type == GB_AT_DNA || ali_type == GB_AT_RNA); // @@@ amino ambiguities not impl
688
689        const uint8_t indices2increment[MAX_TARGET_INDICES+1] = { 0, 0, 6, 4, 3 };
690        for (int i = 0; ambiguous[i]; ++i) {
691            const char *contained = iupac::decode(ambiguous[i], ali_type, false);
692
693            Ambiguity& amb = ambiguity_table[i];
694
695            amb.indices   = strlen(contained);
696            ct_assert(amb.indices>=2 && amb.indices<=MAX_TARGET_INDICES);
697            amb.increment = indices2increment[amb.indices];
698
699            for (int j = 0; j<amb.indices; ++j) {
700                int cidx     = char_to_index_tab[safeCharIndex(contained[j])];
701                ct_assert(cidx<MAX_INDEX_TABLES); // has to be non-ambiguous
702                amb.index[j] = cidx;
703            }
704
705            char_to_index_tab[toupper(ambiguous[i])] =
706            char_to_index_tab[tolower(ambiguous[i])] = i+MAX_INDEX_TABLES;
707        }
708    }
709
710    free(unique_gap_chars);
711    ct_assert(idx==used_bases_tables);
712    initialized = true;
713}
714
715BaseFrequencies::BaseFrequencies(int maxseqlength)
716    : ignore(0)
717{
718    ct_assert(initialized);
719
720    bases_table = new SepBaseFreqPtr[used_bases_tables];
721
722    for (int i=0; i<used_bases_tables; i++) {
723        bases_table[i] = new SepBaseFreq(maxseqlength);
724    }
725
726    sequenceUnits = 0;
727}
728
729void BaseFrequencies::init(int maxseqlength) {
730    int i;
731    for (i=0; i<used_bases_tables; i++) {
732        linear_table(i).init(maxseqlength);
733    }
734
735    sequenceUnits = 0;
736}
737
738void BaseFrequencies::bases_and_gaps_at(int column, int *bases, int *gaps) const {
739    int b = 0;
740    int g = 0;
741
742    for (int i=0; i<used_bases_tables; i++) {
743        char c = upper_index_chars[i];
744
745        if (isGap(c)) {
746            g += table(c)[column];
747        }
748        else {
749            b += table(c)[column];
750        }
751    }
752
753    if (bases) {
754        *bases = b/unitsPerSequence;
755        ct_assert((b%unitsPerSequence) == 0); // could happen if an ambiguous code contains a gap
756    }
757    if (gaps)  {
758        *gaps  = g/unitsPerSequence;
759        ct_assert((g%unitsPerSequence) == 0); // could happen if an ambiguous code contains a gap
760    }
761}
762
763double BaseFrequencies::max_frequency_at(int column, bool ignore_gaps) const {
764    int gaps  = 0;
765    int minus = 0;
766    int maxb  = 0;
767
768    for (int i=0; i<used_bases_tables; i++) {
769        char c = upper_index_chars[i];
770
771        if (isGap(c)) {
772            if (c == '-') {
773                minus += table(c)[column];
774            }
775            else {
776                gaps += table(c)[column];
777            }
778        }
779        else {
780            maxb = std::max(maxb, table(c)[column]);
781        }
782    }
783
784    if (ignore_gaps) {
785        gaps += minus;
786    }
787    else {
788        maxb = std::max(maxb, minus);
789    }
790
791    double mfreq;
792    if (sequenceUnits>gaps) {
793        mfreq = maxb/double(sequenceUnits-gaps);
794    }
795    else { // only gaps (but no '-')
796        ct_assert(sequenceUnits == gaps);
797        mfreq = 0.0;
798    }
799
800#if defined(ASSERTION_USED)
801    if (mfreq != 0.0) {
802        if (ali_type != GB_AT_AA) {
803            ct_assert(mfreq>=0.2); // wrong for amino acids
804        }
805        ct_assert(mfreq<=1.0);
806    }
807#endif
808
809    return mfreq;
810}
811
812BaseFrequencies::~BaseFrequencies() {
813    int i;
814    for (i=0; i<used_bases_tables; i++) {
815        delete bases_table[i];
816    }
817
818    delete [] bases_table;
819}
820
821const PosRange *BaseFrequencies::changed_range(const BaseFrequencies& other) const {
822    int i;
823    int Size = size();
824    int start = Size-1;
825    int end = 0;
826
827    ct_assert(Size==other.size());
828    for (i=0; i<used_bases_tables; i++) {
829        if (linear_table(i).firstDifference(other.linear_table(i), 0, start, &start)) {
830            if (end<start) {
831                end = start;
832            }
833            linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end);
834
835            for (; i<used_bases_tables; i++) {
836                linear_table(i).firstDifference(other.linear_table(i), 0, start-1, &start);
837                if (end<start) {
838                    end = start;
839                }
840                linear_table(i).lastDifference(other.linear_table(i), end+1, Size-1, &end);
841            }
842
843            ct_assert(start<=end);
844
845            static PosRange range;
846            range = PosRange(start, end);
847
848            return &range;
849        }
850    }
851    return NULp;
852}
853void BaseFrequencies::add(const BaseFrequencies& other) {
854    if (other.ignore) return;
855    if (other.size()==0) return;
856    prepare_to_add_elements(other.added_sequences());
857    add(other, 0, other.size()-1);
858}
859void BaseFrequencies::add(const BaseFrequencies& other, int start, int end) {
860    if (other.ignore) return;
861
862    test();
863    other.test();
864
865    ct_assert(start<=end);
866
867    int i;
868    for (i=0; i<used_bases_tables; i++) {
869        linear_table(i).add(other.linear_table(i), start, end);
870    }
871
872    sequenceUnits += other.sequenceUnits;
873
874    test();
875}
876
877void BaseFrequencies::sub(const BaseFrequencies& other) {
878    if (other.ignore) return;
879    sub(other, 0, other.size()-1);
880}
881
882void BaseFrequencies::sub(const BaseFrequencies& other, int start, int end) {
883    if (other.ignore) return;
884
885    test();
886    other.test();
887    ct_assert(start<=end);
888
889    int i;
890    for (i=0; i<used_bases_tables; i++) {
891        linear_table(i).sub(other.linear_table(i), start, end);
892    }
893
894    sequenceUnits -= other.sequenceUnits;
895
896    test();
897}
898
899void BaseFrequencies::sub_and_add(const BaseFrequencies& Sub, const BaseFrequencies& Add, PosRange range) {
900    test();
901    Sub.test();
902    Add.test();
903
904    ct_assert(!Sub.ignore && !Add.ignore);
905    ct_assert(range.is_part());
906
907    int i;
908    for (i=0; i<used_bases_tables; i++) {
909        linear_table(i).sub_and_add(Sub.linear_table(i), Add.linear_table(i), range);
910    }
911
912    test();
913}
914
915const PosRange *BaseFrequencies::changed_range(const char *string1, const char *string2, int min_len) {
916    const unsigned long *range1 = (const unsigned long *)string1;
917    const unsigned long *range2 = (const unsigned long *)string2;
918    const int            step   = sizeof(*range1);
919
920    int l = min_len/step;
921    int r = min_len%step;
922    int i;
923    int j;
924    int k;
925
926    int start = -1, end = -1;
927
928    for (i=0; i<l; i++) {       // search wordwise (it's faster)
929        if (range1[i] != range2[i]) {
930            k = i*step;
931            for (j=0; j<step; j++) {
932                if (string1[k+j] != string2[k+j]) {
933                    start = end = k+j;
934                    break;
935                }
936            }
937
938            break;
939        }
940    }
941
942    k = l*step;
943    if (i==l) {         // no difference found in word -> look at rest
944        for (j=0; j<r; j++) {
945            if (string1[k+j] != string2[k+j]) {
946                start = end = k+j;
947                break;
948            }
949        }
950
951        if (j==r) {     // the strings are equal
952            return NULp;
953        }
954    }
955
956    ct_assert(start != -1 && end != -1);
957
958    for (j=r-1; j>=0; j--) {    // search rest backwards
959        if (string1[k+j] != string2[k+j]) {
960            end = k+j;
961            break;
962        }
963    }
964
965    if (j==-1) {                // not found in rest -> search words backwards
966        int m = start/step;
967        for (i=l-1; i>=m; i--) {
968            if (range1[i] != range2[i]) {
969                k = i*step;
970                for (j=step-1; j>=0; j--) {
971                    if (string1[k+j] != string2[k+j]) {
972                        end = k+j;
973                        break;
974                    }
975                }
976                break;
977            }
978        }
979    }
980
981    ct_assert(start<=end);
982
983    static PosRange range;
984    range = PosRange(start, end);
985
986    return &range;
987}
988
989#define INC_DEC(incdec) do {                                                    \
990        int idx = char_to_index_tab[c];                                         \
991        if (idx<MAX_INDEX_TABLES) {                                             \
992            table(c).incdec(offset, unitsPerSequence);                          \
993        }                                                                       \
994        else {                                                                  \
995            Ambiguity& amb = ambiguity_table[idx-MAX_INDEX_TABLES];             \
996            for (int i = 0; i<amb.indices; ++i) {                               \
997                linear_table(amb.index[i]).incdec(offset, amb.increment);       \
998            }                                                                   \
999        }                                                                       \
1000    } while(0)
1001
1002void BaseFrequencies::incr_short(unsigned char c, int offset) {
1003    INC_DEC(inc_short);
1004}
1005
1006void BaseFrequencies::decr_short(unsigned char c, int offset) {
1007    INC_DEC(dec_short);
1008}
1009
1010void BaseFrequencies::incr_long(unsigned char c, int offset) {
1011    INC_DEC(inc_long);
1012}
1013
1014void BaseFrequencies::decr_long(unsigned char c, int offset) {
1015    INC_DEC(dec_long);
1016}
1017
1018#undef INC_DEC
1019
1020void BaseFrequencies::add(const char *scan_string, int len) {
1021    test();
1022
1023    prepare_to_add_elements(1);
1024
1025    int i;
1026    int sz = size();
1027
1028    if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) {
1029        for (i=0; i<len; i++) {
1030            unsigned char c = scan_string[i];
1031            if (!c) break;
1032            incr_short(c, i);
1033        }
1034        SepBaseFreq& t = table('.');
1035        for (; i<sz; i++) {
1036            t.inc_short(i, unitsPerSequence);
1037        }
1038    }
1039    else {
1040        for (i=0; i<len; i++) {
1041            unsigned char c = scan_string[i];
1042            if (!c) break;
1043            incr_long(c, i);
1044        }
1045        SepBaseFreq& t = table('.');
1046        for (; i<sz; i++) { // LOOP_VECTORIZED[!>=620<7]
1047            t.inc_long(i, unitsPerSequence);
1048        }
1049    }
1050
1051    sequenceUnits += unitsPerSequence;
1052
1053    test();
1054}
1055
1056void BaseFrequencies::sub(const char *scan_string, int len) {
1057    test();
1058
1059    int i;
1060    int sz = size();
1061
1062    if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) {
1063        for (i=0; i<len; i++) {
1064            unsigned char c = scan_string[i];
1065            if (!c) break;
1066            decr_short(c, i);
1067        }
1068        SepBaseFreq& t = table('.');
1069        for (; i<sz; i++) {
1070            t.dec_short(i, unitsPerSequence);
1071        }
1072    }
1073    else {
1074        for (i=0; i<len; i++) {
1075            unsigned char c = scan_string[i];
1076            if (!c) break;
1077            decr_long(c, i);
1078        }
1079        SepBaseFreq& t = table('.');
1080        for (; i<sz; i++) { // LOOP_VECTORIZED[!>=620<7]
1081            t.dec_long(i, unitsPerSequence);
1082        }
1083    }
1084
1085    sequenceUnits -= unitsPerSequence;
1086
1087    test();
1088}
1089
1090void BaseFrequencies::sub_and_add(const char *old_string, const char *new_string, PosRange range) {
1091    test();
1092    int start = range.start();
1093    int end   = range.end();
1094    ct_assert(start<=end);
1095
1096    if (get_table_entry_size()==SHORT_TABLE_ELEM_SIZE) {
1097        for (int i=start; i<=end; i++) {
1098            unsigned char o = old_string[i];
1099            unsigned char n = new_string[i];
1100
1101            ct_assert(o && n); // passed string have to be defined in range
1102
1103            if (o!=n) {
1104                decr_short(o, i);
1105                incr_short(n, i);
1106            }
1107        }
1108    }
1109    else {
1110        for (int i=start; i<=end; i++) {
1111            unsigned char o = old_string[i];
1112            unsigned char n = new_string[i];
1113
1114            ct_assert(o && n); // passed string have to be defined in range
1115
1116            if (o!=n) {
1117                decr_long(o, i);
1118                incr_long(n, i);
1119            }
1120        }
1121    }
1122
1123    test();
1124}
1125
1126void BaseFrequencies::change_table_length(int new_length) {
1127    for (int c=0; c<used_bases_tables; c++) {
1128        linear_table(c).change_table_length(new_length, 0);
1129    }
1130    test();
1131}
1132
1133//  --------------
1134//      tests:
1135
1136#if defined(TEST_CHAR_TABLE_INTEGRITY) || defined(ASSERTION_USED)
1137
1138bool BaseFrequencies::empty() const {
1139    int i;
1140    for (i=0; i<used_bases_tables; i++) {
1141        if (!linear_table(i).empty()) {
1142            return false;
1143        }
1144    }
1145    return true;
1146}
1147
1148bool BaseFrequencies::ok() const {
1149    if (empty()) return true;
1150    if (is_ignored()) return true;
1151
1152    if (sequenceUnits<0) {
1153        fprintf(stderr, "Negative # of sequences (%i) in BaseFrequencies\n", added_sequences());
1154        return false;
1155    }
1156
1157    const int table_size = size();
1158    for (int i=0; i<table_size; i++) {
1159        int bases, gaps;
1160        bases_and_gaps_at(i, &bases, &gaps);
1161
1162        if (!bases && !gaps) {                      // this occurs after insert column
1163            for (int j=i+1; j<table_size; j++) {    // test if we find any bases till end of table !!!
1164
1165                bases_and_gaps_at(j, &bases, NULp);
1166                if (bases) { // bases found -> empty position was illegal
1167                    fprintf(stderr, "Zero in BaseFrequencies at column %i\n", i);
1168                    return false;
1169                }
1170            }
1171            break;
1172        }
1173        else {
1174            if ((bases+gaps)!=added_sequences()) {
1175                fprintf(stderr, "bases+gaps should be equal to # of sequences at column %i (bases=%i, gaps=%i, added_sequences=%i)\n",
1176                        i, bases, gaps, added_sequences());
1177                return false;
1178            }
1179        }
1180    }
1181
1182    return true;
1183}
1184
1185#endif // TEST_CHAR_TABLE_INTEGRITY || ASSERTION_USED
1186
1187#if defined(TEST_CHAR_TABLE_INTEGRITY)
1188
1189void BaseFrequencies::test() const {
1190
1191    if (!ok()) {
1192        GBK_terminate("BaseFrequencies::test() failed");
1193    }
1194}
1195
1196#endif // TEST_CHAR_TABLE_INTEGRITY
1197
1198
1199// --------------------------------------------------------------------------------
1200
1201#ifdef UNIT_TESTS
1202#ifndef TEST_UNIT_H
1203#include <test_unit.h>
1204#endif
1205#include <arbdb.h>
1206
1207#define SETUP(gapChars,alitype) BaseFrequencies::setup(gapChars,alitype)
1208
1209void TEST_char_table() {
1210    for (int useAminoAcids = 0; useAminoAcids<=1; ++useAminoAcids) {
1211        const char *alphabeth      = useAminoAcids ? "ABCDEFGHIJKLMNPQRSTVWXYZ-" : "ACGTN-";
1212        const int   alphabeth_size = strlen(alphabeth);
1213
1214        TEST_EXPECT_EQUAL(alphabeth_size, (useAminoAcids ? 23 : 4)+2); // "+2" is for gap and N/X
1215
1216        const int SEQLEN = 30;
1217        char      seq[SEQLEN+1];
1218
1219        const char *gapChars    = ".-"; // @@@ doesnt make any difference which gaps are used - why?
1220        // const char *gapChars = "-";
1221        // const char *gapChars = ".";
1222        // const char *gapChars = "A"; // makes me fail
1223        SETUP(gapChars, useAminoAcids ? GB_AT_AA : GB_AT_RNA);
1224
1225        if (useAminoAcids) {
1226            // document used amino acid group definition:
1227            TEST_EXPECT_EQUAL(getAAgroupDefinition(), "A,B,C,D,E,F,G,H,I,J,K,L,M,N,P,Q,R,S,T,V,W,Y,Z,X");
1228        }
1229
1230        ConsensusBuildParams BK;
1231        // BK.lower = 70; BK.upper = 95; BK.gapbound = 60; // defaults from awars
1232        BK.lower    = 40; BK.upper = 70; BK.gapbound = 40;
1233
1234        const int LOOPS = 5;
1235        GB_random_seed(100);
1236
1237        int seed[LOOPS];
1238        for (int loop = 0; loop<5; ++loop) {
1239            seed[loop] = GB_random(INT_MAX);
1240        }
1241
1242        for (int loop = 0; loop<5; ++loop) {
1243            GB_random_seed(seed[loop]);
1244
1245            BaseFrequencies tab(SEQLEN);
1246
1247            // build some seq
1248            for (int c = 0; c<SEQLEN; ++c) {
1249                seq[c] = alphabeth[GB_random(alphabeth_size)];
1250            }
1251            seq[SEQLEN] = 0;
1252
1253            TEST_EXPECT_EQUAL(strlen(seq), size_t(SEQLEN));
1254
1255            const int  add_count = 300;
1256            char      *added_seqs[add_count];
1257
1258            const char *s1; // length = SEQLEN
1259            const char *s2;
1260            if (useAminoAcids) {
1261                s1 = "ABCDEFGHIJKLMNPQRSTVWXYZabcdef";
1262                s2 = "ghijklmnpqrstvwxyzABCDEFGHIJKL";
1263            }
1264            else { // rna
1265                s1 = "ACGTACGTAcgtacgtaCGTACGTACGTAC"; // @@@ insert gaps! -> adapt changed results
1266                s2 = "MRWSYKVHDBNmrwsykvhdbnSYKVHDBN"; // use ambiguities
1267            }
1268
1269            // add seq multiple times
1270            for (int a = 0; a<add_count; ++a) {
1271                tab.add(seq, SEQLEN);
1272                added_seqs[a] = strdup(seq);
1273
1274                // modify 1 character in seq:
1275                int sidx  = GB_random(SEQLEN);
1276                int aidx  = GB_random(alphabeth_size);
1277                seq[sidx] = alphabeth[aidx];
1278
1279                if (a == 16) { // with 15 sequences in tab
1280                    // test sub_and_add (short table version)
1281                    char *consensus = tab.build_consensus_string(BK);
1282
1283                    tab.add(s1, SEQLEN);
1284                    PosRange r(0, SEQLEN-1);
1285                    tab.sub_and_add(s1, s2, r);
1286                    tab.sub_and_add(s2, consensus, r);
1287                    tab.sub(consensus, SEQLEN);
1288
1289                    char *consensus2 = tab.build_consensus_string(BK);
1290                    TEST_EXPECT_EQUAL(consensus, consensus2);
1291
1292                    free(consensus2);
1293                    free(consensus);
1294                }
1295            }
1296
1297            // build consensi (just check regression)
1298            // Note: semantic tests for consensus are in ../../NTREE/AP_consensus.cxx@TEST_nucleotide_consensus_and_maxFrequency
1299            {
1300                char *consensus = tab.build_consensus_string(BK);
1301                if (useAminoAcids) {
1302                    switch (seed[loop]) {
1303                        case 1080710223: TEST_EXPECT_EQUAL(consensus, "a.i.dXi.haai.hifdi.d.adXaa..dX"); break;
1304                        case  290469779: TEST_EXPECT_EQUAL(consensus, ".......I.D.A......-f..iXi...dX"); break;
1305                        case  346019116: TEST_EXPECT_EQUAL(consensus, "if.a-hd.hd..a.di-ddf.i..aXIaI."); break;
1306                        case 1708473987: TEST_EXPECT_EQUAL(consensus, "aia..df...haDaadd.d.fha.ii.X.."); break;
1307                        case 2080862653: TEST_EXPECT_EQUAL(consensus, "a.ia.id...a.aah.had......i..d."); break;
1308                        default: TEST_EXPECTATION(all().of(that(consensus).is_equal_to("undef"), that(seed[loop]).is_equal_to(-1))); // default case is unwanted (logs data -> want one case above for each pair)
1309                    }
1310                }
1311                else {
1312                    switch (seed[loop]) {
1313                        case 1080710223: TEST_EXPECT_EQUAL(consensus, "W.gKyc.sR.urSkGcy.kuy......RY-"); break;
1314                        case  290469779: TEST_EXPECT_EQUAL(consensus, ".u-SW.Mgcu.uw-uu..-mmb.sgcRuY."); break;
1315                        case  346019116: TEST_EXPECT_EQUAL(consensus, "y-m.-Suy.c.-u.-g-y.cWgr.cyGrgY"); break;
1316                        case 1708473987: TEST_EXPECT_EQUAL(consensus, ".K.Kc--.cakyw.uW..a-.guuggwH-a"); break;
1317                        case 2080862653: TEST_EXPECT_EQUAL(consensus, ".sgaR..--.uu.cMuuycsk-u.cu.Wc."); break;
1318                        default: TEST_EXPECTATION(all().of(that(consensus).is_equal_to("undef"), that(seed[loop]).is_equal_to(-1))); // default case is unwanted (logs data -> want one case above for each pair)
1319                    }
1320                }
1321
1322                // test sub_and_add()
1323
1324                tab.add(s1, SEQLEN);
1325                PosRange r(0, SEQLEN-1);
1326                tab.sub_and_add(s1, s2, r);
1327                tab.sub_and_add(s2, consensus, r);
1328                tab.sub(consensus, SEQLEN);
1329
1330                char *consensus2 = tab.build_consensus_string(BK);
1331                TEST_EXPECT_EQUAL(consensus, consensus2);
1332
1333                free(consensus2);
1334                free(consensus);
1335            }
1336
1337            for (int a = 0; a<add_count; a++) {
1338                tab.sub(added_seqs[a], SEQLEN);
1339                freenull(added_seqs[a]);
1340            }
1341            {
1342                char *consensus = tab.build_consensus_string(BK);
1343                TEST_EXPECT_EQUAL(consensus, "??????????????????????????????"); // check tab is empty
1344                free(consensus);
1345            }
1346        }
1347    }
1348}
1349
1350void TEST_SepBaseFreq() {
1351    const int LEN  = 20;
1352    const int OFF1 = 6;
1353    const int OFF2 = 7; // adjacent to OFF1
1354
1355    SepBaseFreq short1(LEN), short2(LEN);
1356    for (int i = 0; i<100; ++i) short1.inc_short(OFF1, 1);
1357    for (int i = 0; i<70;  ++i) short1.inc_short(OFF2, 1);
1358    for (int i = 0; i<150; ++i) short2.inc_short(OFF1, 1);
1359    for (int i = 0; i<80;  ++i) short2.inc_short(OFF2, 1);
1360
1361    SepBaseFreq long1(LEN), long2(LEN);
1362    long1.expand_table_entry_size();
1363    long2.expand_table_entry_size();
1364    for (int i = 0; i<2000; ++i) long1.inc_long(OFF1, 1);
1365    for (int i = 0; i<2500; ++i) long2.inc_long(OFF1, 1);
1366
1367    SepBaseFreq shortsum(LEN);
1368    shortsum.add(short1, 0, LEN-1); TEST_EXPECT_EQUAL(shortsum[OFF1], 100); TEST_EXPECT_EQUAL(shortsum[OFF2],  70);
1369    shortsum.add(short2, 0, LEN-1); TEST_EXPECT_EQUAL(shortsum[OFF1], 250); TEST_EXPECT_EQUAL(shortsum[OFF2], 150);
1370    shortsum.sub(short1, 0, LEN-1); TEST_EXPECT_EQUAL(shortsum[OFF1], 150); TEST_EXPECT_EQUAL(shortsum[OFF2],  80);
1371
1372    shortsum.sub_and_add(short1, short2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(shortsum[OFF1], 200); TEST_EXPECT_EQUAL(shortsum[OFF2], 90);
1373
1374    // shortsum.add(long1, 0, LEN-1); // invalid operation -> terminate
1375
1376    SepBaseFreq longsum(LEN);
1377    longsum.expand_table_entry_size();
1378    longsum.add(long1,  0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2000);
1379    longsum.add(long2,  0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 4500);
1380    longsum.sub(long1,  0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2500);
1381    longsum.add(short1, 0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2600);
1382    longsum.sub(short2, 0, LEN-1); TEST_EXPECT_EQUAL(longsum[OFF1], 2450);
1383
1384    longsum.sub_and_add(long1,  long2,  PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 2950);
1385    longsum.sub_and_add(short1, short2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 3000);
1386    longsum.sub_and_add(long1,  short2, PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 1150);
1387    longsum.sub_and_add(short1, long2,  PosRange(0, LEN-1)); TEST_EXPECT_EQUAL(longsum[OFF1], 3550);
1388
1389    int pos = -1;
1390    TEST_EXPECT_EQUAL(short1.firstDifference(short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1);
1391    TEST_EXPECT_EQUAL(short1.lastDifference (short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF2);
1392    TEST_EXPECT_EQUAL(short1.firstDifference(long1,  0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1);
1393    TEST_EXPECT_EQUAL(short1.lastDifference (long1,  0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF2);
1394    TEST_EXPECT_EQUAL(long2.firstDifference (short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1);
1395    TEST_EXPECT_EQUAL(long2.lastDifference  (short2, 0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF2);
1396    TEST_EXPECT_EQUAL(long2.firstDifference (long1,  0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1);
1397    TEST_EXPECT_EQUAL(long2.lastDifference  (long1,  0, LEN-1, &pos), 1); TEST_EXPECT_EQUAL(pos, OFF1);
1398}
1399
1400#endif // UNIT_TESTS
1401
1402// --------------------------------------------------------------------------------
1403
1404
1405
Note: See TracBrowser for help on using the repository browser.