source: branches/alilink/SL/FAST_ALIGNER/fast_aligner.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: 136.4 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : fast_aligner.cxx                                  //
4//   Purpose   : A fast aligner (not a multiple aligner!)          //
5//                                                                 //
6//   Coded by Ralf Westram (coder@reallysoft.de) in 1998           //
7//   Institute of Microbiology (Technical University Munich)       //
8//   http://www.arb-home.de/                                       //
9//                                                                 //
10// =============================================================== //
11
12
13#include "fast_aligner.hxx"
14#include "ClustalV.hxx"
15#include "seq_search.hxx"
16
17#include <island_hopping.h>
18
19#include <awtc_next_neighbours.hxx>
20#include <awt_sel_boxes.hxx>
21
22#include <aw_awars.hxx>
23#include <aw_root.hxx>
24#include <aw_question.hxx>
25
26#include <arbdbt.h>
27#include <ad_cb.h>
28
29#include <arb_defs.h>
30#include <arb_progress.h>
31#include <RangeList.h>
32
33#include <cctype>
34#include <cmath>
35#include <climits>
36
37#include <list>
38#include <awt_config_manager.hxx>
39#include <rootAsWin.h>
40
41// --------------------------------------------------------------------------------
42
43#if defined(DEBUG)
44// #define TRACE_CLUSTAL_DATA
45// #define TRACE_ISLANDHOPPER_DATA
46// #define TRACE_COMPRESSED_ALIGNMENT
47// #define TRACE_RELATIVES
48#endif // DEBUG
49
50// --------------------------------------------------------------------------------
51
52enum FA_report {
53    FA_NO_REPORT,               // no report
54    FA_TEMP_REPORT,             // report to temporary entries
55    FA_REPORT,                  // report to resident entries
56};
57
58enum FA_range {
59    FA_WHOLE_SEQUENCE,          // align whole sequence
60    FA_AROUND_CURSOR,           // align xxx positions around current cursor position
61    FA_SELECTED_RANGE,          // align selected range
62    FA_SAI_MULTI_RANGE,         // align versus multi range define by SAI
63};
64
65enum FA_turn {
66    FA_TURN_NEVER,              // never try to turn sequence
67    FA_TURN_INTERACTIVE,        // try to turn, but query user
68    FA_TURN_ALWAYS,             // turn if score is better
69};
70
71enum FA_reference {
72    FA_REF_EXPLICIT,            // reference sequence explicitly specified
73    FA_REF_CONSENSUS,           // use group consensus as reference
74    FA_REF_RELATIVES,           // search next relatives by PT server
75};
76
77enum FA_alignTarget {
78    FA_CURRENT,                 // align current species
79    FA_MARKED,                  // align all marked species
80    FA_SELECTED,                // align selected species (= range)
81};
82
83enum FA_errorAction {
84    FA_NO_ACTION,                // do nothing
85    FA_MARK_FAILED,              // mark failed species (unmark rest)
86    FA_MARK_ALIGNED,             // mark aligned species (unmark rest)
87};
88
89struct AlignParams {
90    FA_report report;
91    bool      showGapsMessages; // display messages about missing gaps in master?
92    PosRange  range;            // range to be aligned
93};
94
95struct SearchRelativeParams : virtual Noncopyable {
96    FamilyFinder *ff;
97    char         *pt_server_alignment; // alignment used in pt_server (may differ from 'alignment')
98    int           maxRelatives;        // max # of relatives to use
99
100    SearchRelativeParams(FamilyFinder *ff_, const char *pt_server_alignment_, int maxRelatives_)
101        : ff(ff_),
102          pt_server_alignment(strdup(pt_server_alignment_)),
103          maxRelatives(maxRelatives_)
104    {}
105   
106    ~SearchRelativeParams() {
107        free(pt_server_alignment);
108        delete(ff);
109    }
110
111    FamilyFinder *getFamilyFinder() { return ff; }
112};
113
114// --------------------------------------------------------------------------------
115
116#define GAP_CHAR     '-'
117#define QUALITY_NAME "ASC_ALIGNER_CLIENT_SCORE"
118#define INSERTS_NAME "AMI_ALIGNER_MASTER_INSERTS"
119
120#define FA_AWAR_ROOT                "faligner/"
121#define FA_AWAR_TO_ALIGN            FA_AWAR_ROOT "what"
122#define FA_AWAR_REFERENCE           FA_AWAR_ROOT "against"
123#define FA_AWAR_REFERENCE_NAME      FA_AWAR_ROOT "sagainst"
124#define FA_AWAR_RANGE               FA_AWAR_ROOT "range"
125#define FA_AWAR_PROTECTION          FA_AWAR_ROOT "protection"
126#define FA_AWAR_AROUND              FA_AWAR_ROOT "around"
127#define FA_AWAR_MIRROR              FA_AWAR_ROOT "mirror"
128#define FA_AWAR_REPORT              FA_AWAR_ROOT "report"
129#define FA_AWAR_SHOW_GAPS_MESSAGES  FA_AWAR_ROOT "show_gaps"
130#define FA_AWAR_CONTINUE_ON_ERROR   FA_AWAR_ROOT "continue_on_error"
131#define FA_AWAR_ACTION_ON_ERROR     FA_AWAR_ROOT "action_on_error"
132#define FA_AWAR_USE_SECONDARY       FA_AWAR_ROOT "use_secondary"
133#define FA_AWAR_NEXT_RELATIVES      FA_AWAR_ROOT "next_relatives"
134#define FA_AWAR_RELATIVE_RANGE      FA_AWAR_ROOT "relrange"
135#define FA_AWAR_PT_SERVER_ALIGNMENT "tmp/" FA_AWAR_ROOT "relative_ali"
136#define FA_AWAR_SAI_RANGE_NAME      FA_AWAR_ROOT "sai/sainame"
137#define FA_AWAR_SAI_RANGE_CHARS     FA_AWAR_ROOT "sai/chars"
138
139#define FA_AWAR_ISLAND_HOPPING_ROOT  "island_hopping/"
140#define FA_AWAR_USE_ISLAND_HOPPING   FA_AWAR_ISLAND_HOPPING_ROOT "use"
141#define FA_AWAR_ESTIMATE_BASE_FREQ   FA_AWAR_ISLAND_HOPPING_ROOT "estimate_base_freq"
142#define FA_AWAR_BASE_FREQ_A          FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_a"
143#define FA_AWAR_BASE_FREQ_C          FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_c"
144#define FA_AWAR_BASE_FREQ_G          FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_g"
145#define FA_AWAR_BASE_FREQ_T          FA_AWAR_ISLAND_HOPPING_ROOT "base_freq_t"
146#define FA_AWAR_SUBST_PARA_AC        FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ac"
147#define FA_AWAR_SUBST_PARA_AG        FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ag"
148#define FA_AWAR_SUBST_PARA_AT        FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_at"
149#define FA_AWAR_SUBST_PARA_CG        FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_cg"
150#define FA_AWAR_SUBST_PARA_CT        FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_ct"
151#define FA_AWAR_SUBST_PARA_GT        FA_AWAR_ISLAND_HOPPING_ROOT "subst_para_gt"
152#define FA_AWAR_EXPECTED_DISTANCE    FA_AWAR_ISLAND_HOPPING_ROOT "expected_dist"
153#define FA_AWAR_STRUCTURE_SUPPLEMENT FA_AWAR_ISLAND_HOPPING_ROOT "struct_suppl"
154#define FA_AWAR_THRESHOLD            FA_AWAR_ISLAND_HOPPING_ROOT "threshold"
155#define FA_AWAR_GAP_A                FA_AWAR_ISLAND_HOPPING_ROOT "gapa"
156#define FA_AWAR_GAP_B                FA_AWAR_ISLAND_HOPPING_ROOT "gapb"
157#define FA_AWAR_GAP_C                FA_AWAR_ISLAND_HOPPING_ROOT "gapc"
158
159// --------------------------------------------------------------------------------
160
161static IslandHopping *island_hopper = NULp; // NULp -> use fast aligner; else use island hopper
162
163static GB_alignment_type global_alignmentType = GB_AT_UNKNOWN;  // type of actually aligned sequence
164
165static int currentSequenceNumber;    // used for counter
166static int overallSequenceNumber;
167
168// --------------------------------------------------------------------------------
169
170inline ARB_ERROR species_not_found(GB_CSTR species_name) {
171    return GBS_global_string("No species '%s' found!", species_name);
172}
173
174static ARB_ERROR reverseComplement(GBDATA *gb_species, GB_CSTR ali, int max_protection) {
175    GBDATA    *gbd   = GBT_find_sequence(gb_species, ali);
176    ARB_ERROR  error = NULp;
177
178    if (!gbd) {
179        error = GBS_global_string("No 'data' found for species '%s'", GBT_read_name(gb_species));
180    }
181    else {
182        int my_protection = GB_read_security_write(gbd);
183
184        if (my_protection<=max_protection) { // ok
185            char              *seq      = GB_read_string(gbd);
186            int                length   = GB_read_string_count(gbd);
187            GBDATA            *gb_main  = GB_get_root(gb_species);
188            GB_alignment_type  ali_type = GBT_get_alignment_type(gb_main, ali);
189
190            char T_or_U;
191            error = GBT_determine_T_or_U(ali_type, &T_or_U, "reverse-complement");
192            if (!error) {
193                GBT_reverseComplementNucSequence(seq, length, T_or_U);
194                error = GB_write_string(gbd, seq);
195            }
196        }
197        else { // protection error
198            error = GBS_global_string("Cannot reverse-complement species '%s' because of protection level", GBT_read_name(gb_species));
199        }
200
201    }
202
203    return error;
204}
205
206static void build_reverse_complement(AW_window *aw, const AlignDataAccess *data_access) {
207    GBDATA *gb_main = data_access->gb_main;
208
209    GB_push_transaction(gb_main);
210
211    AW_root        *root              = aw->get_root();
212    FA_alignTarget  revComplWhat      = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
213    char           *default_alignment = GBT_get_default_alignment(gb_main);
214    GB_CSTR         alignment         = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
215    ARB_ERROR       error             = NULp;
216    int             max_protection    = root->awar(FA_AWAR_PROTECTION)->read_int();
217
218    switch (revComplWhat) {
219        case FA_CURRENT: { // current species
220            GB_CSTR species_name = root->awar(AWAR_SPECIES_NAME)->read_string();
221            GBDATA *gb_species = GBT_find_species(gb_main, species_name);
222            if (!gb_species) error = species_not_found(species_name);
223            if (!error) error = reverseComplement(gb_species, alignment, max_protection);
224            break;
225        }
226        case FA_MARKED: { // marked species
227            GBDATA *gb_species = GBT_first_marked_species(gb_main);
228
229            if (!gb_species) {
230                error = "There is no marked species";
231            }
232            while (gb_species) {
233                error = reverseComplement(gb_species, alignment, max_protection);
234                if (error) break;
235                gb_species = GBT_next_marked_species(gb_species);
236            }
237            break;
238        }
239        case FA_SELECTED: { // selected species (editor selection!)
240           
241            Aligner_get_first_selected_species get_first_selected_species = data_access->get_first_selected_species;
242            Aligner_get_next_selected_species  get_next_selected_species  = data_access->get_next_selected_species;
243
244            int     count      = 0;
245            GBDATA *gb_species = get_first_selected_species(&count);
246           
247            if (!gb_species) {
248                error = "There is no selected species!";
249            }
250            while (gb_species) {
251                error = reverseComplement(gb_species, alignment, max_protection);
252                if (error) break;
253                gb_species = get_next_selected_species();
254            }
255            break;
256        }
257        default: {
258            fa_assert(0);
259            break;
260        }
261    }
262
263    GB_end_transaction_show_error(gb_main, error, aw_message);
264}
265
266// --------------------------------------------------------------------------------
267
268class AliChange { // describes a local alignment change
269    const CompactedSubSequence& Old;
270    const CompactedSubSequence& New;
271
272public:
273    AliChange(const CompactedSubSequence& old_, const CompactedSubSequence& new_)
274        : Old(old_), New(new_)
275    {
276        fa_assert(Old.may_refer_to_same_part_as(New));
277    }
278
279    int follow(ExplicitRange& range) const {
280        ExplicitRange compressed(Old.compPosition(range.start()),
281                                 Old.compPosition(range.end()));
282
283        int exp_start = New.expdPosition(compressed.start());
284        int exp_end   = New.expdPosition(compressed.end());
285
286        int gaps_before = New.no_of_gaps_before(compressed.start());
287        int gaps_after = New.no_of_gaps_after(compressed.end());
288
289        fa_assert(gaps_before >= 0);
290        fa_assert(gaps_after >= 0);
291
292        range = ExplicitRange(exp_start-gaps_before, exp_end+gaps_after);
293
294        return compressed.size(); // number of bases
295    }
296};
297
298class LooseBases {
299    typedef std::list<ExplicitRange> Ranges;
300
301    Ranges ranges;
302   
303public:
304
305    bool is_empty() const { return ranges.empty(); }
306    void clear() { ranges.clear(); }
307
308    void memorize(ExplicitRange range) {
309        ranges.push_front(range);
310    }
311    ExplicitRange recall() {
312        ExplicitRange range = ranges.front();
313        ranges.pop_front();
314        return range;
315    }
316    int follow_ali_change(const AliChange& change) {
317        // transform positions according to changed alignment (oldSequence -> newSequence)
318        // returns the number of bases contained in 'this'
319        int basecount = 0;
320        if (!is_empty()) {
321            for (Ranges::iterator r = ranges.begin(); r != ranges.end(); ++r) {
322                basecount += change.follow(*r);
323            }
324        }
325        return basecount;
326    }
327    void append(LooseBases& loose) { ranges.splice(ranges.end(), loose.ranges); }
328    int follow_ali_change_and_append(LooseBases& loose, const AliChange& change) {
329        // returns the number of loose bases (added from 'loose')
330        int basecount = loose.follow_ali_change(change);
331        append(loose);
332        return basecount;
333    }
334};
335
336static LooseBases unaligned_bases; // if fast_align cannot align (no master bases) -> stores positions here
337
338
339static const char *read_name(GBDATA *gbd) {
340    return gbd ? GBT_read_name(gbd) : "GROUP-CONSENSUS";
341}
342
343inline int relatedBases(char base1, char base2) {
344    return baseMatch(base1, base2)==1;
345}
346
347inline char alignQuality(char slave, char master) {
348    fa_assert(slave);
349    fa_assert(master);
350    char result = '#';
351    if (slave==master)              result = '-';   // equal
352    else if (slave==GAP_CHAR)           result = '+';   // inserted gap
353    else if (master==GAP_CHAR)          result = '+';   // no gap in master
354    else if (relatedBases(slave, master))   result = '~';   // mutation (related bases)
355    return result;                          // mutation (non-related bases)
356}
357
358// -------------------------
359//      Debugging stuff
360
361#ifdef DEBUG
362static char *lstr(const char *s, int len) {
363    static int alloc;
364    static char *buffer;
365
366    if (alloc<(len+1)) {
367        if (alloc) free(buffer);
368        ARB_alloc(buffer, alloc=len+100);
369    }
370
371    memcpy(buffer, s, len);
372    buffer[len] = 0;
373
374    return buffer;
375}
376
377#define BUFLEN 120
378
379inline char compareChar(char base1, char base2) {
380    return base1==base2 ? '=' : (relatedBases(base1, base2) ? 'x' : 'X');
381}
382
383#if defined(TRACE_COMPRESSED_ALIGNMENT)
384   
385static void dump_n_compare_one(const char *seq1, const char *seq2, long len, long offset) {
386    fa_assert(len<=BUFLEN);
387    char compare[BUFLEN+1];
388
389    for (long l=0; l<len; l++) {
390        compare[l] = (is_ali_gap(seq1[l]) && is_ali_gap(seq2[l])) ? ' ' : compareChar(seq1[l], seq2[l]);
391    }
392
393    compare[len] = 0;
394
395    printf(" %li '%s'\n", offset, lstr(seq1, len));
396    printf(" %li '%s'\n", offset, lstr(seq2, len));
397    printf(" %li '%s'\n", offset, compare);
398}
399
400inline void dump_rest(const char *seq, long len, int idx, long offset) {
401    printf(" Rest von Sequenz %i:\n", idx);
402    while (len>BUFLEN) {
403        printf(" %li '%s'\n", offset, lstr(seq, BUFLEN));
404        seq += BUFLEN;
405        len -= BUFLEN;
406        offset += BUFLEN;
407    }
408
409    fa_assert(len>0);
410    printf(" '%s'\n", lstr(seq, len));
411}
412
413static void dump_n_compare(const char *text, const char *seq1, long len1, const char *seq2, long len2) {
414    long offset = 0;
415
416    printf(" Comparing %s:\n", text);
417
418    while (len1>0 && len2>0) {
419        long done = 0;
420
421        if (len1>=BUFLEN && len2>=BUFLEN) {
422            dump_n_compare_one(seq1, seq2, done=BUFLEN, offset);
423        }
424        else {
425            long min = len1<len2 ? len1 : len2;
426            dump_n_compare_one(seq1, seq2, done=min, offset);
427        }
428
429        seq1 += done;
430        seq2 += done;
431        len1 -= done;
432        len2 -= done;
433        offset += done;
434    }
435
436    if (len1>0) dump_rest(seq1, len1, 1, offset);
437    if (len2>0) dump_rest(seq2, len2, 2, offset);
438    printf(" -------------------\n");
439}
440
441static void dump_n_compare(const char *text, const CompactedSubSequence& seq1, const CompactedSubSequence& seq2) {
442    dump_n_compare(text, seq1.text(), seq1.length(), seq2.text(), seq2.length());
443}
444#endif // TRACE_COMPRESSED_ALIGNMENT
445
446#undef BUFLEN
447
448inline void dumpSeq(const char *seq, long len, long pos) {
449    printf("'%s' ", lstr(seq, len));
450    printf("(Pos=%li,Len=%li)", pos, len);
451}
452
453#define dump()                                                          \
454    do {                                                                \
455        double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength); \
456                                                                        \
457        printf(" Score = %li (Significance=%f)\n"                       \
458               " Master = ", bestScore, sig);                           \
459        dumpSeq(bestMasterLeft.text(), bestLength, bestMasterLeft.leftOf()); \
460        printf("\n"                                                     \
461               " Slave  = ");                                           \
462        dumpSeq(bestSlaveLeft.text(), bestLength, bestSlaveLeft.leftOf()); \
463        printf("\n");                                                   \
464    } while (0)
465
466#endif //DEBUG
467
468
469// ------------------------------------------------
470//      INLINE-functions used in fast_align():
471
472inline double log3(double d) {
473    return log(d)/log(3.0);
474}
475inline double partSignificance(long seq1len, long seq2len, long partlen) {
476    // returns log3 of significance of the part
477    // usage: partSignificance(...) < log3(maxAllowedSignificance)
478    return log3((seq1len-partlen)*(seq2len-partlen)) - partlen;
479}
480
481inline ARB_ERROR bufferTooSmall() {
482    return "Cannot align - reserved buffer is to small";
483}
484
485inline long insertsToNextBase(AlignBuffer *alignBuffer, const SequencePosition& master) {
486    int inserts;
487    int nextBase;
488
489    if (master.rightOf()>0) {
490        nextBase = master.expdPosition();
491    }
492    else {
493        nextBase = master.sequence().expdPosition(master.sequence().length());
494    }
495    inserts = nextBase-alignBuffer->offset();
496
497    return inserts;
498}
499
500inline void insertBase(AlignBuffer *alignBuffer,
501                       SequencePosition& master, SequencePosition& slave,
502                       FastAlignReport *report)
503{
504    char slaveBase = *slave.text();
505    char masterBase = *master.text();
506
507    alignBuffer->set(slaveBase, alignQuality(slaveBase, masterBase));
508    report->count_aligned_base(slaveBase!=masterBase);
509    ++slave;
510    ++master;
511}
512
513inline void insertSlaveBases(AlignBuffer *alignBuffer,
514                             SequencePosition& slave,
515                             int length,
516                             FastAlignReport *report)
517{
518    alignBuffer->copy(slave.text(), alignQuality(*slave.text(), GAP_CHAR), length);
519    report->count_unaligned_base(length);
520    slave += length;
521}
522
523inline void insertGap(AlignBuffer *alignBuffer,
524                      SequencePosition& master,
525                      FastAlignReport *report)
526{
527    char masterBase = *master.text();
528
529    alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, masterBase));
530    report->count_aligned_base(masterBase!=GAP_CHAR);
531    ++master;
532}
533
534static ARB_ERROR insertClustalValigned(AlignBuffer *alignBuffer,
535                                       SequencePosition& master,
536                                       SequencePosition& slave,
537                                       const char *masterAlignment, const char *slaveAlignment, long alignmentLength,
538                                       FastAlignReport *report)
539{
540    // inserts bases of 'slave' to 'alignBuffer' according to alignment in 'masterAlignment' and 'slaveAlignment'
541#define ACID '*'    // contents of 'masterAlignment' and 'slaveAlignment'
542#define GAP  '-'
543
544    int pos;
545    int baseAtLeft = 0;     // TRUE -> last position in alignBuffer contains a base
546
547    for (pos=0; pos<alignmentLength; pos++) {
548        long insert = insertsToNextBase(alignBuffer, master); // in expanded seq
549
550        if (masterAlignment[pos]==ACID) {
551            if (insert>0) {
552                if (insert>alignBuffer->free()) return bufferTooSmall();
553                alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
554                fa_assert(insertsToNextBase(alignBuffer, master)==0);
555                insert = 0;
556            }
557
558            if (!alignBuffer->free()) return bufferTooSmall();
559            if (slaveAlignment[pos]==ACID) {
560                insertBase(alignBuffer, master, slave, report);
561                baseAtLeft = 1;
562            }
563            else {
564                insertGap(alignBuffer, master, report);
565                baseAtLeft = 0;
566            }
567        }
568        else {
569            int slave_bases;
570
571            fa_assert(masterAlignment[pos]==GAP);
572            for (slave_bases=1; pos+slave_bases<alignmentLength && masterAlignment[pos+slave_bases]==GAP; slave_bases++) {
573                ; // count gaps in master (= # of slave bases to insert)
574            }
575            if (!baseAtLeft && insert>slave_bases) {
576                int ins_gaps = insert-slave_bases;
577
578                fa_assert(alignBuffer->free()>=ins_gaps);
579                alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), ins_gaps);
580                insert -= ins_gaps;
581            }
582
583            if (insert<slave_bases) { // master has less gaps than slave bases to insert
584                report->memorize_insertion(master.expdPosition(), slave_bases-insert);
585            }
586            else if (insert>slave_bases) { // master has more gaps than slave bases to insert
587                fa_assert(baseAtLeft);
588            }
589
590            unaligned_bases.memorize(ExplicitRange(slave.expdPosition(), // memorize base positions without counterpart in master
591                                                   slave.expdPosition(slave_bases-1)));
592
593            if (slave_bases>alignBuffer->free()) return bufferTooSmall();
594            insertSlaveBases(alignBuffer, slave, slave_bases, report);
595            pos += slave_bases-1; // -1 compensates for()-increment
596            baseAtLeft = 1;
597        }
598    }
599
600    return NULp;
601
602#undef GAP
603#undef ACID
604}
605
606static ARB_ERROR insertAligned(AlignBuffer *alignBuffer,
607                               SequencePosition& master, SequencePosition& slave, long partLength,
608                               FastAlignReport *report)
609{
610    // insert bases of 'slave' to 'alignBuffer' according to 'master'
611    if (partLength) {
612        long insert = insertsToNextBase(alignBuffer, master);
613
614        fa_assert(partLength>=0);
615
616        if (insert<0) { // insert gaps into master
617            long min_insert = insert;
618
619            report->memorize_insertion(master.expdPosition(), -insert);
620
621            while (insert<0 && partLength) {
622                if (insert<min_insert) min_insert = insert;
623                if (!alignBuffer->free()) {
624                    return bufferTooSmall();
625                }
626                insertBase(alignBuffer, master, slave, report);
627                partLength--;
628                insert = insertsToNextBase(alignBuffer, master);
629            }
630
631            fa_assert(partLength>=0);
632            if (partLength==0) { // all inserted
633                return NULp;
634            }
635        }
636
637        fa_assert(insert>=0);
638
639        if (insert>0) { // insert gaps into slave
640            if (insert>alignBuffer->free()) return bufferTooSmall();
641            alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
642            fa_assert(insertsToNextBase(alignBuffer, master)==0);
643        }
644
645        fa_assert(partLength>=0);
646
647        while (partLength--) {
648            insert = insertsToNextBase(alignBuffer, master);
649
650            fa_assert(insert>=0);
651            if (insert>0) {
652                if (insert>=alignBuffer->free()) return bufferTooSmall();
653                alignBuffer->set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), insert);
654            }
655            else {
656                if (!alignBuffer->free()) {
657                    return bufferTooSmall();
658                }
659            }
660
661            insertBase(alignBuffer, master, slave, report);
662        }
663    }
664
665    return NULp;
666}
667
668static ARB_ERROR cannot_fast_align(const CompactedSubSequence& master, long moffset, long mlength,
669                                   const CompactedSubSequence& slaveSequence, long soffset, long slength,
670                                   int max_seq_length,
671                                   AlignBuffer *alignBuffer,
672                                   FastAlignReport *report)
673{
674    const char *mtext = master.text(moffset);
675    const char *stext = slaveSequence.text(soffset);
676    ARB_ERROR   error = NULp;
677
678    if (slength) {
679        if (mlength) { // if slave- and master-sequences contain bases, we call the slow aligner
680#ifdef TRACE_CLUSTAL_DATA
681            printf("ClustalV-Align:\n");
682            printf(" mseq = '%s'\n", lstr(mtext, mlength));
683            printf(" sseq = '%s'\n", lstr(stext, slength));
684#endif // TRACE_CLUSTAL_DATA
685
686            int is_dna = -1;
687
688            switch (global_alignmentType) {
689                case GB_AT_RNA:
690                case GB_AT_DNA: is_dna = 1; break;
691                case GB_AT_AA:  is_dna = 0; break;
692                default: error = "Unknown alignment type - aligner aborted"; break;
693            }
694
695            const char *maligned, *saligned;
696            int         len;
697            if (!error) {
698                int score; // unused
699                error = ClustalV_align(is_dna,
700                                       1,
701                                       mtext, mlength, stext, slength,
702                                       master.gapsBefore(moffset), 
703                                       max_seq_length,
704                                       maligned, saligned, len, score);
705            }
706
707            if (!error) {
708#ifdef TRACE_CLUSTAL_DATA
709                printf("ClustalV returns:\n");
710                printf(" maligned = '%s'\n", lstr(maligned, len));
711                printf(" saligned = '%s'\n", lstr(saligned, len));
712#endif // TRACE_CLUSTAL_DATA
713
714                SequencePosition masterPos(master, moffset);
715                SequencePosition slavePos(slaveSequence, soffset);
716
717                error = insertClustalValigned(alignBuffer, masterPos, slavePos, maligned, saligned, len, report);
718
719#if (defined(DEBUG) && 0)
720
721                SequencePosition master2(master->sequence(), moffset);
722                SequencePosition slave2(slaveSequence, soffset);
723                char *cmp = new char[len];
724
725                for (int l=0; l<len; l++) {
726                    int gaps = 0;
727
728                    if (maligned[l]=='*') {
729                        maligned[l] = *master2.text();
730                        ++master2;
731                    }
732                    else {
733                        gaps++;
734                    }
735
736                    if (saligned[l]=='*') {
737                        saligned[l] = *slave2.text();
738                        ++slave2;
739                    }
740                    else {
741                        gaps++;
742                    }
743
744                    cmp[l] = gaps || maligned[l]==saligned[l] ? '=' : 'X';
745                }
746
747                printf(" master = '%s'\n", lstr(maligned, len));
748                printf(" slave  = '%s'\n", lstr(saligned, len));
749                printf("          '%s'\n", lstr(cmp, len));
750
751                delete [] cmp;
752#endif
753            }
754        }
755        else { // master is empty here, so we just copy in the slave bases
756            if (slength<=alignBuffer->free()) {
757                unaligned_bases.memorize(ExplicitRange(slaveSequence.expdPosition(soffset),
758                                                       slaveSequence.expdPosition(soffset+slength-1)));
759               
760                alignBuffer->copy(slaveSequence.text(soffset), '?', slength);   // '?' means not aligned (just inserted)
761                // corrected by at alignBuffer->correctUnalignedPositions()
762                report->count_unaligned_base(slength);
763            }
764            else {
765                error = bufferTooSmall();
766            }
767        }
768    }
769
770    return error;
771}
772
773// ------------------------------------
774//      #define's for fast_align()
775
776#define TEST_BETTER_SCORE()                                             \
777    do {                                                                \
778        if (score>bestScore) {                                          \
779            bestScore = score;                                          \
780            bestLength = masterRight.text() - masterLeft.text();        \
781            bestMasterLeft = masterLeft;                                \
782            bestSlaveLeft = slaveLeft;                                  \
783        }                                                               \
784    } while (0)
785
786#define CAN_SCORE_LEFT()    (masterLeft.leftOf() && slaveLeft.leftOf())
787#define CAN_SCORE_RIGHT()   (masterRight.rightOf() && slaveRight.rightOf())
788
789#define SCORE_LEFT()                                                    \
790    do {                                                                \
791        score += *(--masterLeft).text()==*(--slaveLeft).text() ? match : mismatch; \
792        TEST_BETTER_SCORE();                                            \
793    } while (0)
794
795#define SCORE_RIGHT()                                                   \
796    do {                                                                \
797        score += *(++masterRight).text()==*(++slaveRight).text() ? match : mismatch; \
798        TEST_BETTER_SCORE();                                            \
799    } while (0)
800
801
802ARB_ERROR FastSearchSequence::fast_align(const CompactedSubSequence& slaveSequence,
803                                         AlignBuffer *alignBuffer,
804                                         int max_seq_length,
805                                         int match, int mismatch,
806                                         FastAlignReport *report) const
807{
808    // aligns 'slaveSequence' to 'this'
809    //
810    // returns
811    //     NULp     -> all ok -> 'alignBuffer' contains aligned sequence
812    //     or error -> failure -> no results
813
814    ARB_ERROR error   = NULp;
815    int       aligned = 0;
816
817    // set the following #if to zero to test ClustalV-Aligner (not fast_aligner)
818#if 1
819
820    static double lowSignificance;
821    static int lowSignificanceInitialized;
822
823    if (!lowSignificanceInitialized) {
824        lowSignificance = log3(0.01);
825        lowSignificanceInitialized = 1;
826    }
827
828    SequencePosition slave(slaveSequence);
829    long bestScore=0;
830    SequencePosition bestMasterLeft(sequence());
831    SequencePosition bestSlaveLeft(slaveSequence);
832    long bestLength=0;
833
834    while (slave.rightOf()>=3) { // with all triples of slaveSequence
835        FastSearchOccurrence occurrence(*this, slave.text()); // "search" first
836        SequencePosition rightmostSlave = slave;
837
838        while (occurrence.found()) { // with all occurrences of the triple
839            long score = match*3;
840            SequencePosition masterLeft(occurrence.sequence(), occurrence.offset());
841            SequencePosition masterRight(occurrence.sequence(), occurrence.offset()+3);
842            SequencePosition slaveLeft(slave);
843            SequencePosition slaveRight(slave, 3);
844
845            while (score>0) {
846                if (CAN_SCORE_LEFT()) {
847                    SCORE_LEFT();
848                }
849                else {
850                    while (score>0 && CAN_SCORE_RIGHT()) {
851                        SCORE_RIGHT();
852                    }
853                    break;
854                }
855
856                if (CAN_SCORE_RIGHT()) {
857                    SCORE_RIGHT();
858                }
859                else {
860                    while (score>0 && CAN_SCORE_LEFT()) {
861                        SCORE_LEFT();
862                    }
863                    break;
864                }
865            }
866
867            occurrence.gotoNext(); // "search" next
868
869            if (rightmostSlave<slaveRight) {
870                rightmostSlave = slaveRight;
871                rightmostSlave -= 3;
872            }
873        }
874
875        if (rightmostSlave>slave)   slave = rightmostSlave;
876        else                ++slave;
877    }
878
879    if (bestLength) {
880        double sig = partSignificance(sequence().length(), slaveSequence.length(), bestLength);
881
882        if (sig<lowSignificance) {
883            long masterLeftOf = bestMasterLeft.leftOf();
884            long masterRightStart = masterLeftOf+bestLength;
885            long masterRightOf = bestMasterLeft.rightOf()-bestLength;
886            long slaveLeftOf = bestSlaveLeft.leftOf();
887            long slaveRightStart = slaveLeftOf+bestLength;
888            long slaveRightOf = bestSlaveLeft.rightOf()-bestLength;
889
890#define MIN_ALIGNMENT_RANGE 4
891
892            if (!error) {
893                if (masterLeftOf >= MIN_ALIGNMENT_RANGE && slaveLeftOf >= MIN_ALIGNMENT_RANGE) {
894                    CompactedSubSequence   leftCompactedMaster(sequence(), 0, masterLeftOf);
895                    FastSearchSequence     leftMaster(leftCompactedMaster);
896
897                    error = leftMaster.fast_align(CompactedSubSequence(slaveSequence, 0, slaveLeftOf),
898                                                  alignBuffer, max_seq_length, match, mismatch, report);
899                }
900                else if (slaveLeftOf>0) {
901                    error = cannot_fast_align(sequence(), 0, masterLeftOf,
902                                              slaveSequence, 0, slaveLeftOf,
903                                              max_seq_length, alignBuffer, report);
904                }
905
906                aligned = 1;
907            }
908
909            // align part of slave sequence according to master sequence:
910
911            if (!error) {
912#if (defined(DEBUG) && 0)
913                long offset = alignBuffer->offset();
914                long used;
915#endif
916                error = insertAligned(alignBuffer, bestMasterLeft, bestSlaveLeft, bestLength, report);
917#if (defined(DEBUG) && 0)
918                used = alignBuffer->offset()-offset;
919                printf("aligned '%s' (len=%li, address=%li)\n", lstr(alignBuffer->text()+offset, used), used, long(alignBuffer));
920#endif
921                aligned = 1;
922            }
923
924            if (!error) {
925                if (masterRightOf >= MIN_ALIGNMENT_RANGE && slaveRightOf >= MIN_ALIGNMENT_RANGE) {
926                    CompactedSubSequence rightCompactedMaster(sequence(), masterRightStart, masterRightOf);
927                    FastSearchSequence rightMaster(rightCompactedMaster);
928
929                    error = rightMaster.fast_align(CompactedSubSequence(slaveSequence, slaveRightStart, slaveRightOf),
930                                                   alignBuffer,
931                                                   max_seq_length, match, mismatch, report);
932                }
933                else if (slaveRightOf>0) {
934                    error = cannot_fast_align(sequence(), masterRightStart, masterRightOf,
935                                              slaveSequence, slaveRightStart, slaveRightOf,
936                                              max_seq_length, alignBuffer, report);
937                }
938
939                aligned = 1;
940            }
941
942        }
943    }
944
945#endif
946
947    if (!aligned && !error) {
948        error = cannot_fast_align(sequence(), 0, sequence().length(),
949                                  slaveSequence, 0, slaveSequence.length(),
950                                  max_seq_length, alignBuffer, report);
951    }
952
953    return error;
954}
955
956#undef dump
957#undef TEST_BETTER_SCORE
958#undef CAN_SCORE_LEFT
959#undef CAN_SCORE_RIGHT
960#undef SCORE_LEFT
961#undef SCORE_RIGHT
962
963inline long calcSequenceChecksum(const char *data, long length) {
964    return GB_checksum(data, length, 1, GAP_CHARS);
965}
966
967#if defined(WARN_TODO)
968#warning firstColumn + lastColumn -> PosRange
969#endif
970
971static CompactedSubSequence *readCompactedSequence(GBDATA      *gb_species,
972                                                   const char  *ali,
973                                                   ARB_ERROR   *errorPtr,
974                                                   char       **dataPtr,     // if dataPtr is specified, it will be set to the WHOLE uncompacted sequence
975                                                   long        *seqChksum,   // may be NULp (of part of sequence)
976                                                   PosRange     range)       // if !range.is_whole() -> return only part of the sequence
977{
978    ARB_ERROR             error = NULp;
979    GBDATA               *gbd;
980    CompactedSubSequence *seq   = NULp;
981
982    gbd = GBT_find_sequence(gb_species, ali);       // get sequence
983
984    if (gbd) {
985        long length = GB_read_string_count(gbd);
986        char *data = GB_read_string(gbd);
987        long partLength;
988        char *partData;
989
990        if (dataPtr) {                  // make a copy of the whole uncompacted sequence (returned to caller)
991            *dataPtr = data;
992        }
993
994        int firstColumn = range.start();
995        if (range.is_part()) {     // take only part of sequence
996            int lastColumn = range.end();
997
998            fa_assert(firstColumn>=0);
999            fa_assert(firstColumn<=lastColumn);
1000
1001            // include all surrounding gaps (@@@ this might cause problems with partial alignment)
1002            while (firstColumn>0 && is_ali_gap(data[firstColumn-1])) {
1003                firstColumn--;
1004            }
1005            if (lastColumn!=-1) {
1006                while (lastColumn<(length-1) && is_ali_gap(data[lastColumn+1])) lastColumn++;
1007                fa_assert(lastColumn<length);
1008            }
1009
1010            partData = data+firstColumn;
1011            int slen = length-firstColumn;
1012
1013            fa_assert(slen>=0);
1014            fa_assert((size_t)slen==strlen(partData));
1015
1016            if (lastColumn==-1) { // take all till end of sequence
1017                partLength = slen;
1018            }
1019            else {
1020                partLength = lastColumn-firstColumn+1;
1021                if (partLength>slen) partLength = slen;     // cut rest, if we have any
1022            }
1023        }
1024        else {
1025            partLength = length;
1026            partData = data;
1027        }
1028
1029        if (!error) {
1030            if (seqChksum) {
1031                *seqChksum = calcSequenceChecksum(partData, partLength);
1032            }
1033
1034            seq = new CompactedSubSequence(partData, partLength, GBT_read_name(gb_species), firstColumn);
1035        }
1036
1037        if (!dataPtr) {     // free sequence only if user has not requested to get it
1038            free(data);
1039        }
1040    }
1041    else {
1042        error = GBS_global_string("No 'data' found for species '%s'", GBT_read_name(gb_species));
1043        if (dataPtr) *dataPtr = NULp; // (user must not care to free data if we fail)
1044    }
1045
1046    fa_assert(errorPtr);
1047    *errorPtr = error;
1048
1049    return seq;
1050}
1051
1052static ARB_ERROR writeStringToAlignment(GBDATA *gb_species, GB_CSTR alignment, GB_CSTR data_name, GB_CSTR str, bool temporary) {
1053    GBDATA    *gb_ali  = GB_search(gb_species, alignment, GB_DB);
1054    ARB_ERROR  error   = NULp;
1055    GBDATA    *gb_name = GB_search(gb_ali, data_name, GB_STRING);
1056
1057    if (gb_name) {
1058        fa_assert(GB_is_ancestor_of(gb_ali, gb_name));
1059        error = GB_write_string(gb_name, str);
1060        if (temporary && !error) error = GB_set_temporary(gb_name);
1061    }
1062    else {
1063        error = GBS_global_string("Cannot create entry '%s' for '%s'", data_name, GBT_read_name(gb_species));
1064    }
1065
1066    return error;
1067}
1068
1069// --------------------------------------------------------------------------------
1070
1071static ARB_ERROR alignCompactedTo(CompactedSubSequence     *toAlignSequence,
1072                                  const FastSearchSequence *alignTo,
1073                                  int                       max_seq_length,
1074                                  GB_CSTR                   alignment,
1075                                  long                      toAlignChksum,
1076                                  GBDATA                   *gb_toAlign,
1077                                  GBDATA                   *gb_alignTo, // may be NULp
1078                                  const AlignParams&        ali_params)
1079{
1080    // if only part of the sequence should be aligned, then this functions already gets only the part
1081    // (i.o.w.: toAlignSequence, alignTo and toAlignChksum refer to the partial sequence)
1082    AlignBuffer alignBuffer(max_seq_length);
1083    if (ali_params.range.start()>0) {
1084        alignBuffer.set(GAP_CHAR, alignQuality(GAP_CHAR, GAP_CHAR), ali_params.range.start());
1085    }
1086
1087    const char *master_name = read_name(gb_alignTo);
1088
1089    FastAlignReport report(master_name, ali_params.showGapsMessages);
1090
1091    {
1092        static GBDATA *last_gb_toAlign = NULp;
1093        if (gb_toAlign!=last_gb_toAlign) {
1094            last_gb_toAlign = gb_toAlign;
1095            currentSequenceNumber++;
1096        }
1097    }
1098
1099#ifdef TRACE_COMPRESSED_ALIGNMENT
1100    printf("alignCompactedTo(): master='%s' ", master_name);
1101    printf("slave='%s'\n", toAlignSequence->name());
1102#endif // TRACE_COMPRESSED_ALIGNMENT
1103
1104    ARB_ERROR error = GB_pop_transaction(gb_toAlign);
1105    if (!error) {
1106        if (island_hopper) {
1107            error = island_hopper->do_align();
1108            if (!error) {
1109                fa_assert(island_hopper->was_aligned());
1110
1111#ifdef TRACE_ISLANDHOPPER_DATA
1112                printf("Island-Hopper returns:\n");
1113                printf("maligned = '%s'\n", lstr(island_hopper->get_result_ref(), island_hopper->get_result_length()));
1114                printf("saligned = '%s'\n", lstr(island_hopper->get_result(), island_hopper->get_result_length()));
1115#endif // TRACE_ISLANDHOPPER_DATA
1116
1117                SequencePosition masterPos(alignTo->sequence(), 0);
1118                SequencePosition slavePos(*toAlignSequence, 0);
1119
1120                error = insertClustalValigned(&alignBuffer, masterPos, slavePos, island_hopper->get_result_ref(), island_hopper->get_result(), island_hopper->get_result_length(), &report);
1121            }
1122        }
1123        else {
1124            error = alignTo->fast_align(*toAlignSequence, &alignBuffer, max_seq_length, 2, -10, &report); // <- here we align
1125        }
1126    }
1127
1128    if (!error) {
1129        alignBuffer.correctUnalignedPositions();
1130        if (alignBuffer.free()) {
1131            alignBuffer.set('-', alignQuality(GAP_CHAR, GAP_CHAR), alignBuffer.free()); // rest of alignBuffer is set to '-'
1132        }
1133        alignBuffer.restoreDots(*toAlignSequence);
1134    }
1135
1136#ifdef TRACE_COMPRESSED_ALIGNMENT
1137    if (!error) {
1138        CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1139        dump_n_compare("reference vs. aligned:", alignTo->sequence(), alignedSlave);
1140    }
1141#endif // TRACE_COMPRESSED_ALIGNMENT
1142
1143    {
1144        GB_ERROR err = GB_push_transaction(gb_toAlign);
1145        if (!error) error = err;
1146    }
1147
1148    if (!error) {
1149        if (calcSequenceChecksum(alignBuffer.text(), alignBuffer.length())!=toAlignChksum) { // sequence-chksum changed
1150            error = "Internal aligner error (sequence checksum changed) -- aborted";
1151
1152#ifdef TRACE_COMPRESSED_ALIGNMENT
1153            CompactedSubSequence alignedSlave(alignBuffer.text(), alignBuffer.length(), toAlignSequence->name());
1154            dump_n_compare("Old Slave vs. new Slave", *toAlignSequence, alignedSlave);
1155#endif // TRACE_COMPRESSED_ALIGNMENT
1156        }
1157        else {
1158            {
1159                GB_topSecurityLevel  unsecured(gb_toAlign);
1160                GBDATA              *gbd = GBT_add_data(gb_toAlign, alignment, "data", GB_STRING);
1161
1162                if (!gbd) {
1163                    error = "Can't find/create sequence data";
1164                }
1165                else {
1166                    if (ali_params.range.is_part()) { // we aligned just a part of the sequence
1167                        char *buffer       = GB_read_string(gbd); // so we read old sequence data
1168                        long  len          = GB_read_string_count(gbd);
1169                        if (!buffer) error = GB_await_error();
1170                        else {
1171                            int  lenToCopy   = ali_params.range.size();
1172                            long wholeChksum = calcSequenceChecksum(buffer, len);
1173
1174                            memcpy(buffer+ali_params.range.start(), alignBuffer.text()+ali_params.range.start(), lenToCopy); // copy in the aligned part
1175                            // @@@ genau um 1 position zu niedrig
1176
1177                            if (calcSequenceChecksum(buffer, len) != wholeChksum) {
1178                                error            = "Internal aligner error (sequence checksum changed) -- aborted";
1179# ifdef TRACE_COMPRESSED_ALIGNMENT
1180                                char *buffer_org = GB_read_string(gbd);
1181                                dump_n_compare("Old seq vs. new seq (slave)", buffer_org, len, buffer, len);
1182                                free(buffer_org);
1183# endif // TRACE_COMPRESSED_ALIGNMENT
1184                            }
1185                            else {
1186                                GB_write_string(gbd, "");
1187                                error = GB_write_string(gbd, buffer);
1188                            }
1189                        }
1190
1191                        free(buffer);
1192                    }
1193                    else {
1194                        alignBuffer.setDotsAtEOSequence();
1195                        error = GB_write_string(gbd, alignBuffer.text()); // aligned all -> write all
1196                    }
1197                }
1198            }
1199
1200            if (!error && ali_params.report != FA_NO_REPORT) {
1201                // create temp-entry for slave containing alignment quality:
1202
1203                error = writeStringToAlignment(gb_toAlign, alignment, QUALITY_NAME, alignBuffer.quality(), ali_params.report==FA_TEMP_REPORT);
1204                if (!error) {
1205                    // create temp-entry for master containing insert dots:
1206
1207                    int   buflen    = max_seq_length*2;
1208                    char *buffer    = ARB_alloc<char>(buflen+1);
1209                    char *afterLast = buffer;
1210
1211                    if (!buffer) {
1212                        error = "out of memory";
1213                    }
1214                    else {
1215                        memset(buffer, '-', buflen);
1216                        buffer[buflen] = 0;
1217
1218                        const FastAlignInsertion *inserts = report.insertion();
1219                        while (inserts) {
1220                            memset(buffer+inserts->offset(), '>', inserts->gaps());
1221                            afterLast = buffer+inserts->offset()+inserts->gaps();
1222                            inserts = inserts->next();
1223                        }
1224
1225                        *afterLast = 0;
1226                        if (gb_alignTo) {
1227                            error = writeStringToAlignment(gb_alignTo, alignment, INSERTS_NAME, buffer, ali_params.report==FA_TEMP_REPORT);
1228                        }
1229                    }
1230                }
1231            }
1232        }
1233    }
1234    return error;
1235}
1236
1237ARB_ERROR FastAligner_delete_temp_entries(GBDATA *gb_species, const char *alignment) {
1238    fa_assert(gb_species);
1239    fa_assert(alignment);
1240
1241    GBDATA    *gb_ali = GB_search(gb_species, alignment, GB_FIND);
1242    ARB_ERROR  error  = NULp;
1243
1244    if (gb_ali) {
1245        GBDATA *gb_name = GB_search(gb_ali, QUALITY_NAME, GB_FIND);
1246        if (gb_name) {
1247            error = GB_delete(gb_name);
1248#if defined(DEBUG)
1249            printf("deleted '%s' of '%s' (gb_name=%p)\n", QUALITY_NAME, GBT_read_name(gb_species), gb_name);
1250#endif
1251        }
1252
1253        if (!error) {
1254            gb_name = GB_search(gb_ali, INSERTS_NAME, GB_FIND);
1255            if (gb_name) {
1256                error = GB_delete(gb_name);
1257#if defined(DEBUG)
1258                printf("deleted '%s' of '%s' (gb_name=%p)\n", INSERTS_NAME, GBT_read_name(gb_species), gb_name);
1259#endif
1260            }
1261        }
1262    }
1263
1264    return error;
1265}
1266
1267static ARB_ERROR align_error(ARB_ERROR olderr, GBDATA *gb_toAlign, GBDATA *gb_alignTo) {
1268    // used by alignTo() and alignToNextRelative() to transform errors coming from subroutines
1269    // can be used by higher functions
1270
1271    const char *name_toAlign = read_name(gb_toAlign);
1272    const char *name_alignTo = read_name(gb_alignTo);
1273
1274    fa_assert(olderr);
1275
1276    return GBS_global_string("Error while aligning '%s' to '%s':\n%s",
1277                             name_toAlign, name_alignTo, olderr.deliver());
1278}
1279
1280static ARB_ERROR alignTo(GBDATA                   *gb_toAlign,
1281                         GB_CSTR                   alignment,
1282                         const FastSearchSequence *alignTo,
1283                         GBDATA                   *gb_alignTo, // may be NULp
1284                         int                       max_seq_length,
1285                         const AlignParams&        ali_params)
1286{
1287    ARB_ERROR error  = NULp;
1288    long      chksum = -1;
1289
1290    CompactedSubSequence *toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULp, &chksum, ali_params.range);
1291
1292    if (island_hopper) {
1293        GBDATA *gb_seq = GBT_find_sequence(gb_toAlign, alignment);      // get sequence
1294        if (gb_seq) {
1295            long        length = GB_read_string_count(gb_seq);
1296            const char *data   = GB_read_char_pntr(gb_seq);
1297
1298            island_hopper->set_toAlign_sequence(data);
1299            island_hopper->set_alignment_length(length);
1300        }
1301    }
1302
1303
1304
1305    if (!error) {
1306        error = alignCompactedTo(toAlignSequence, alignTo, max_seq_length, alignment, chksum, gb_toAlign, gb_alignTo, ali_params);
1307        if (error) error = align_error(error, gb_toAlign, gb_alignTo);
1308        delete toAlignSequence;
1309    }
1310
1311    return error;
1312}
1313
1314static ARB_ERROR alignToGroupConsensus(GBDATA                     *gb_toAlign,
1315                                       GB_CSTR                     alignment,
1316                                       Aligner_get_consensus_func  get_consensus,
1317                                       int                         max_seq_length,
1318                                       const AlignParams&          ali_params)
1319{
1320    ARB_ERROR  error     = NULp;
1321    char      *consensus = get_consensus(read_name(gb_toAlign), ali_params.range);
1322    size_t     cons_len  = strlen(consensus);
1323
1324    fa_assert(cons_len);
1325
1326    for (size_t i = 0; i<cons_len; ++i) { // translate consensus to be accepted by aligner
1327        switch (consensus[i]) {
1328            case '=': consensus[i] = '-'; break;
1329            default: break;
1330        }
1331    }
1332
1333    CompactedSubSequence compacted(consensus, cons_len, "group consensus", ali_params.range.start());
1334
1335    {
1336        FastSearchSequence fast(compacted);
1337        error = alignTo(gb_toAlign, alignment, &fast, NULp, max_seq_length, ali_params);
1338    }
1339
1340    free(consensus);
1341
1342    return error;
1343}
1344
1345static void appendNameAndUsedBasePositions(char **toString, GBDATA *gb_species, int usedBasePositions) {
1346    // if usedBasePositions == -1 -> unknown;
1347
1348    char *currInfo;
1349    if (usedBasePositions<0) {
1350        currInfo = strdup(GBT_read_name(gb_species));
1351    }
1352    else {
1353        fa_assert(usedBasePositions>0); // otherwise it should NOT be announced here!
1354        currInfo = GBS_global_string_copy("%s:%i", GBT_read_name(gb_species), usedBasePositions);
1355    }
1356
1357    char *newString = NULp;
1358    if (*toString) {
1359        newString = GBS_global_string_copy("%s, %s", *toString, currInfo);
1360    }
1361    else {
1362        newString = currInfo;
1363        currInfo  = NULp; // don't free
1364    }
1365
1366    freeset(*toString, newString);
1367    free(currInfo);
1368}
1369
1370inline int min(int i, int j) { return i<j ? i : j; }
1371
1372static ARB_ERROR alignToNextRelative(SearchRelativeParams&  relSearch,
1373                                     int                    max_seq_length,
1374                                     FA_turn                turnAllowed,
1375                                     GB_CSTR                alignment,
1376                                     GBDATA                *gb_toAlign,
1377                                     const AlignParams&     ali_params)
1378{
1379    CompactedSubSequence *toAlignSequence = NULp;
1380    bool use_different_pt_server_alignment = 0 != strcmp(relSearch.pt_server_alignment, alignment);
1381
1382    ARB_ERROR   error;
1383    long        chksum;
1384    int         relativesToTest = relSearch.maxRelatives*2; // get more relatives from pt-server (needed when use_different_pt_server_alignment == true)
1385    char      **nearestRelative = new char*[relativesToTest+1];
1386    int         next_relatives;
1387    int         i;
1388    GBDATA     *gb_main         = GB_get_root(gb_toAlign);
1389
1390    if (use_different_pt_server_alignment) {
1391        turnAllowed = FA_TURN_NEVER; // makes no sense if we're using a different alignment for the pt_server
1392    }
1393
1394    for (next_relatives=0; next_relatives<relativesToTest; next_relatives++) {
1395        nearestRelative[next_relatives] = NULp;
1396    }
1397    next_relatives = 0;
1398
1399    bool restart = true;
1400    while (restart) {
1401        restart = false;
1402
1403        char *findRelsBySeq = NULp;
1404        if (use_different_pt_server_alignment) {
1405            toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULp, &chksum, ali_params.range);
1406
1407            GBDATA *gbd = GBT_find_sequence(gb_toAlign, relSearch.pt_server_alignment); // use a different alignment for next relative search
1408            if (!gbd) {
1409                error = GBS_global_string("Species '%s' has no data in alignment '%s'", GBT_read_name(gb_toAlign), relSearch.pt_server_alignment);
1410            }
1411            else {
1412                findRelsBySeq = GB_read_string(gbd);
1413            }
1414        }
1415        else {
1416            toAlignSequence = readCompactedSequence(gb_toAlign, alignment, &error, &findRelsBySeq, &chksum, ali_params.range);
1417        }
1418
1419        if (error) {
1420            delete toAlignSequence;
1421            return error; // @@@ leaks ?
1422        }
1423
1424        while (next_relatives) {
1425            next_relatives--;
1426            freenull(nearestRelative[next_relatives]);
1427        }
1428
1429        {
1430            // find relatives
1431            FamilyFinder    *familyFinder = relSearch.getFamilyFinder();
1432            const PosRange&  range        = familyFinder->get_TargetRange();
1433
1434            if (range.is_part()) {
1435                range.copy_corresponding_part(findRelsBySeq, findRelsBySeq, strlen(findRelsBySeq));
1436                turnAllowed = FA_TURN_NEVER; // makes no sense if we're using partial relative search
1437            }
1438
1439            error = familyFinder->searchFamily(findRelsBySeq, FF_FORWARD, relativesToTest+1, 0); // @@@ make min_score configurable
1440
1441            // @@@ case where no relative found (due to min score) handle how ? abort ? warn ?
1442           
1443            double bestScore = 0;
1444            if (!error) {
1445#if defined(DEBUG)
1446                double lastScore = -1;
1447#if defined(TRACE_RELATIVES)
1448                fprintf(stderr, "List of relatives used for '%s':\n", GBT_read_name(gb_toAlign));
1449#endif // TRACE_RELATIVES
1450#endif // DEBUG
1451                for (const FamilyList *fl = familyFinder->getFamilyList(); fl; fl=fl->next) {
1452                    if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1453                        if (GBT_find_species(gb_main, fl->name)) { // @@@
1454                            double thisScore = familyFinder->uses_rel_matches() ? fl->rel_matches : fl->matches;
1455#if defined(DEBUG)
1456                            // check whether family list is sorted correctly
1457                            fa_assert(lastScore < 0 || lastScore >= thisScore);
1458                            lastScore        = thisScore;
1459#if defined(TRACE_RELATIVES)
1460                            fprintf(stderr, "- %s (%5.2f)\n", fl->name, thisScore);
1461#endif // TRACE_RELATIVES
1462#endif // DEBUG
1463                            if (thisScore>=bestScore) bestScore = thisScore;
1464                            if (next_relatives<(relativesToTest+1)) {
1465                                nearestRelative[next_relatives] = strdup(fl->name);
1466                                next_relatives++;
1467                            }
1468                        }
1469                    }
1470                }
1471            }
1472
1473            if (!error && turnAllowed != FA_TURN_NEVER) { // test if mirrored sequence has better relatives
1474                char   *mirroredSequence  = strdup(findRelsBySeq);
1475                long    length            = strlen(mirroredSequence);
1476                double  bestMirroredScore = 0;
1477
1478                char T_or_U;
1479                error = GBT_determine_T_or_U(global_alignmentType, &T_or_U, "reverse-complement");
1480                if (!error) {
1481                    GBT_reverseComplementNucSequence(mirroredSequence, length, T_or_U);
1482                    error = familyFinder->searchFamily(mirroredSequence, FF_FORWARD, relativesToTest+1, 0); // @@@ make min_score configurable
1483                }
1484                if (!error) {
1485#if defined(DEBUG)
1486                    double lastScore = -1;
1487#if defined(TRACE_RELATIVES)
1488                    fprintf(stderr, "List of relatives used for '%s' (turned seq):\n", GBT_read_name(gb_toAlign));
1489#endif // TRACE_RELATIVES
1490#endif // DEBUG
1491                    for (const FamilyList *fl = familyFinder->getFamilyList(); fl; fl = fl->next) {
1492                        double thisScore = familyFinder->uses_rel_matches() ? fl->rel_matches : fl->matches;
1493#if defined(DEBUG)
1494                        // check whether family list is sorted correctly
1495                        fa_assert(lastScore < 0 || lastScore >= thisScore);
1496                        lastScore = thisScore;
1497#if defined(TRACE_RELATIVES)
1498                        fprintf(stderr, "- %s (%5.2f)\n", fl->name, thisScore);
1499#endif // TRACE_RELATIVES
1500#endif // DEBUG
1501                        if (thisScore >= bestMirroredScore) {
1502                            if (strcmp(toAlignSequence->name(), fl->name)!=0) {
1503                                if (GBT_find_species(gb_main, fl->name)) bestMirroredScore = thisScore; // @@@
1504                            }
1505                        }
1506                    }
1507                }
1508
1509                int turnIt = 0;
1510                if (bestMirroredScore>bestScore) {
1511                    if (turnAllowed==FA_TURN_INTERACTIVE) {
1512                        const char *message;
1513                        if (familyFinder->uses_rel_matches()) {
1514                            message = GBS_global_string("'%s' seems to be the other way round (score: %.1f%%, score if turned: %.1f%%)",
1515                                                        toAlignSequence->name(), bestScore*100, bestMirroredScore*100);
1516                        }
1517                        else {
1518                            message = GBS_global_string("'%s' seems to be the other way round (score: %li, score if turned: %li)",
1519                                                        toAlignSequence->name(), long(bestScore+.5), long(bestMirroredScore+.5));
1520                        }
1521                        turnIt = aw_question("fastali_turn_sequence", message, "Turn sequence,Leave sequence alone")==0;
1522                    }
1523                    else {
1524                        fa_assert(turnAllowed == FA_TURN_ALWAYS);
1525                        turnIt = 1;
1526#if defined(TRACE_RELATIVES)
1527                        fprintf(stderr, "Using turned sequence!\n");
1528#endif // TRACE_RELATIVES
1529                    }
1530                }
1531
1532                if (turnIt) { // write mirrored sequence
1533                    GBDATA *gbd = GBT_find_sequence(gb_toAlign, alignment);
1534                    GB_topSecurityLevel unsecured(gbd);
1535                    error = GB_write_string(gbd, mirroredSequence);
1536                    if (!error) {
1537                        delete toAlignSequence;
1538                        restart = true; // continue while loop
1539                    }
1540                }
1541
1542                free(mirroredSequence);
1543            }
1544        }
1545        free(findRelsBySeq);
1546    }
1547
1548    if (!error) {
1549        if (!next_relatives) {
1550            char warning[200];
1551            sprintf(warning, "No relative found for '%s'", toAlignSequence->name());
1552            aw_message(warning);
1553        }
1554        else {
1555            // assuming relatives are sorted! (nearest to farthest)
1556
1557            // get data pointers
1558            typedef GBDATA *GBDATAP;
1559            GBDATAP *gb_reference = new GBDATAP[relSearch.maxRelatives];
1560
1561            {
1562                for (i=0; i<relSearch.maxRelatives && i<next_relatives; i++) {
1563                    GBDATA *gb_species = GBT_find_species(gb_main, nearestRelative[i]);
1564                    if (!gb_species) { // pt-server seems not up to date!
1565                        error = species_not_found(nearestRelative[i]);
1566                        break;
1567                    }
1568
1569                    GBDATA *gb_sequence = GBT_find_sequence(gb_species, alignment);
1570                    if (gb_sequence) { // if relative has sequence data in the current alignment ..
1571                        gb_reference[i] = gb_species;
1572                    }
1573                    else { // remove this relative
1574                        free(nearestRelative[i]);
1575                        for (int j = i+1; j<next_relatives; ++j) {
1576                            nearestRelative[j-1] = nearestRelative[j];
1577                        }
1578                        next_relatives--;
1579                        nearestRelative[next_relatives] = NULp;
1580                        i--; // redo same index
1581                    }
1582                }
1583
1584                // delete superfluous relatives
1585                for (; i<next_relatives; ++i) freenull(nearestRelative[i]);
1586
1587                if (next_relatives>relSearch.maxRelatives) {
1588                    next_relatives = relSearch.maxRelatives;
1589                }
1590            }
1591
1592            // align
1593
1594            if (!error) {
1595                CompactedSubSequence *alignToSequence = readCompactedSequence(gb_reference[0], alignment, &error, NULp, NULp, ali_params.range);
1596                fa_assert(alignToSequence);
1597
1598                if (island_hopper) {
1599                    GBDATA *gb_ref   = GBT_find_sequence(gb_reference[0], alignment); // get reference sequence
1600                    GBDATA *gb_align = GBT_find_sequence(gb_toAlign, alignment);      // get sequence to align
1601
1602                    if (gb_ref && gb_align) {
1603                        long        length_ref   = GB_read_string_count(gb_ref);
1604                        const char *data;
1605
1606                        data = GB_read_char_pntr(gb_ref);
1607                        island_hopper->set_ref_sequence(data);
1608
1609                        data = GB_read_char_pntr(gb_align);
1610                        island_hopper->set_toAlign_sequence(data);
1611
1612                        island_hopper->set_alignment_length(length_ref);
1613                    }
1614                }
1615
1616                {
1617                    FastSearchSequence referenceFastSeq(*alignToSequence);
1618
1619                    error = alignCompactedTo(toAlignSequence, &referenceFastSeq,
1620                                             max_seq_length, alignment, chksum,
1621                                             gb_toAlign, gb_reference[0], ali_params);
1622                }
1623
1624                if (error) {
1625                    error = align_error(error, gb_toAlign, gb_reference[0]);
1626                }
1627                else {
1628                    char *used_relatives = NULp;
1629
1630                    if (unaligned_bases.is_empty()) {
1631                        appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], -1);
1632                    }
1633                    else {
1634                        if (island_hopper) {
1635                            appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], -1);
1636                            if (next_relatives>1) error = "Island hopping uses only one relative";
1637                        }
1638                        else {
1639                            LooseBases loose;
1640                            LooseBases loose_for_next_relative;
1641
1642                            int unaligned_positions;
1643                            {
1644                                CompactedSubSequence *alignedSequence = readCompactedSequence(gb_toAlign, alignment, &error, NULp, NULp, ali_params.range);
1645
1646                                unaligned_positions = loose.follow_ali_change_and_append(unaligned_bases, AliChange(*toAlignSequence, *alignedSequence));
1647                                // now loose holds the unaligned (and recalculated) parts from last relative
1648                                delete alignedSequence;
1649                            }
1650
1651                            if (!error) {
1652                                int toalign_positions = toAlignSequence->length();
1653                                if (unaligned_positions<toalign_positions) {
1654                                    appendNameAndUsedBasePositions(&used_relatives, gb_reference[0], toalign_positions-unaligned_positions);
1655                                }
1656                            }
1657
1658                            for (i=1; i<next_relatives && !error; i++) {
1659                                loose.append(loose_for_next_relative);
1660                                int unaligned_positions_for_next = 0;
1661
1662                                while (!loose.is_empty() && !error) {
1663                                    ExplicitRange         partRange(intersection(loose.recall(), ali_params.range));
1664                                    CompactedSubSequence *alignToPart = readCompactedSequence(gb_reference[i], alignment, &error, NULp, NULp, partRange);
1665
1666                                    if (!error) {
1667                                        long                  part_chksum;
1668                                        CompactedSubSequence *toAlignPart = readCompactedSequence(gb_toAlign, alignment, &error, NULp, &part_chksum, partRange);
1669
1670                                        fa_assert(contradicted(error, toAlignPart));
1671
1672                                        if (!error) {
1673                                            AlignParams loose_ali_params = { ali_params.report, ali_params.showGapsMessages, partRange };
1674
1675                                            FastSearchSequence referenceFastSeq(*alignToPart);
1676                                            error = alignCompactedTo(toAlignPart, &referenceFastSeq,
1677                                                                     max_seq_length, alignment, part_chksum,
1678                                                                     gb_toAlign, gb_reference[i], loose_ali_params);
1679                                            if (!error) {
1680                                                CompactedSubSequence *alignedPart = readCompactedSequence(gb_toAlign, alignment, &error, NULp, NULp, partRange);
1681                                                if (!error) {
1682                                                    unaligned_positions_for_next += loose_for_next_relative.follow_ali_change_and_append(unaligned_bases, AliChange(*toAlignPart, *alignedPart));
1683                                                }
1684                                                delete alignedPart;
1685                                            }
1686                                        }
1687                                        delete toAlignPart;
1688                                    }
1689                                    delete alignToPart;
1690                                }
1691
1692                                if (!error) {
1693                                    fa_assert(unaligned_positions_for_next <= unaligned_positions); // means: number of unaligned positions has increased by use of relative
1694                                    if (unaligned_positions_for_next<unaligned_positions) {
1695                                        appendNameAndUsedBasePositions(&used_relatives, gb_reference[i], unaligned_positions-unaligned_positions_for_next);
1696                                        unaligned_positions = unaligned_positions_for_next;
1697                                    }
1698                                }
1699                            }
1700                        }
1701                    }
1702
1703                    if (!error) { // write used relatives to db-field
1704                        error = GBT_write_string(gb_toAlign, "used_rels", used_relatives);
1705                    }
1706                    free(used_relatives);
1707                }
1708
1709                delete alignToSequence;
1710            }
1711
1712            delete [] gb_reference;
1713        }
1714    }
1715
1716    delete toAlignSequence;
1717
1718    for (i=0; i<next_relatives; i++) freenull(nearestRelative[i]);
1719    delete [] nearestRelative;
1720
1721    return error;
1722}
1723
1724// ------------------------
1725//      AlignmentReference
1726
1727class AlignmentReference : virtual Noncopyable {
1728    GB_CSTR            alignment;
1729    int                max_seq_length;
1730    const AlignParams& ali_params;
1731
1732public:
1733    AlignmentReference(GB_CSTR            alignment_,
1734                       int                max_seq_length_,
1735                       const AlignParams& ali_params_)
1736        : alignment(alignment_),
1737          max_seq_length(max_seq_length_),
1738          ali_params(ali_params_)
1739    {}
1740    virtual ~AlignmentReference() {}
1741
1742    virtual ARB_ERROR align_to(GBDATA *gb_toalign) const = 0;
1743
1744    GB_CSTR get_alignment() const { return alignment; }
1745    int get_max_seq_length() const { return max_seq_length; }
1746    const AlignParams& get_ali_params() const { return ali_params; }
1747};
1748
1749
1750#if defined(WARN_TODO)
1751#warning make alignTo a member of ExplicitReference (or of AlignmentReference)
1752#warning let alignToGroupConsensus and alignToNextRelative use ExplicitReference
1753#endif
1754
1755class ExplicitReference: public AlignmentReference { // derived from a Noncopyable
1756    const FastSearchSequence *targetSequence;
1757    GBDATA                   *gb_alignTo;
1758
1759public:
1760    ExplicitReference(GB_CSTR                   alignment_,
1761                      const FastSearchSequence *targetSequence_,
1762                      GBDATA                   *gb_alignTo_, 
1763                      int                       max_seq_length_,
1764                      const AlignParams&        ali_params_)
1765        : AlignmentReference(alignment_, max_seq_length_, ali_params_), 
1766          targetSequence(targetSequence_), 
1767          gb_alignTo(gb_alignTo_) 
1768    {}
1769
1770    ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1771        return alignTo(gb_toalign, get_alignment(), targetSequence, gb_alignTo, get_max_seq_length(), get_ali_params());
1772    }
1773};
1774
1775#if defined(WARN_TODO)
1776#warning make alignToGroupConsensus a member of ConsensusReference
1777#endif
1778
1779class ConsensusReference: public AlignmentReference {
1780    Aligner_get_consensus_func  get_consensus;
1781
1782public:
1783    ConsensusReference(GB_CSTR                     alignment_, 
1784                       Aligner_get_consensus_func  get_consensus_, 
1785                       int                         max_seq_length_, 
1786                       const AlignParams&          ali_params_)
1787        : AlignmentReference(alignment_, max_seq_length_, ali_params_),
1788          get_consensus(get_consensus_) 
1789    {}
1790   
1791    ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1792        return alignToGroupConsensus(gb_toalign, get_alignment(), get_consensus, get_max_seq_length(), get_ali_params());
1793    }
1794};
1795
1796#if defined(WARN_TODO)
1797#warning make alignToNextRelative a member of SearchRelativesReference
1798#endif
1799
1800class SearchRelativesReference: public AlignmentReference {
1801    SearchRelativeParams&  relSearch;
1802    FA_turn                turnAllowed;
1803
1804public:
1805    SearchRelativesReference(SearchRelativeParams&  relSearch_,
1806                             int                    max_seq_length_,
1807                             FA_turn                turnAllowed_,
1808                             GB_CSTR                alignment_,
1809                             const AlignParams&     ali_params_)
1810        : AlignmentReference(alignment_, max_seq_length_, ali_params_),
1811          relSearch(relSearch_),
1812          turnAllowed(turnAllowed_) 
1813    {}
1814
1815    ARB_ERROR align_to(GBDATA *gb_toalign) const OVERRIDE {
1816        return alignToNextRelative(relSearch, get_max_seq_length(), turnAllowed, get_alignment(), gb_toalign, get_ali_params());
1817    }
1818};
1819
1820
1821// ----------------
1822//      Aligner
1823
1824class Aligner : virtual Noncopyable {
1825    GBDATA *gb_main;
1826
1827    // define alignment target(s):
1828    FA_alignTarget                     alignWhat;
1829    GB_CSTR                            alignment;                  // name of alignment to use (if NULp -> use default)
1830    GB_CSTR                            toalign;                    // name of species to align (used if alignWhat == FA_CURRENT)
1831    Aligner_get_first_selected_species get_first_selected_species; // used if alignWhat == FA_SELECTED
1832    Aligner_get_next_selected_species  get_next_selected_species;  // --- "" ---
1833
1834    // define reference sequence(s):
1835    GB_CSTR                    reference;     // name of reference species
1836    Aligner_get_consensus_func get_consensus; // function to get consensus seq
1837    SearchRelativeParams&      relSearch;     // params to search for relatives
1838
1839    // general params:
1840    FA_turn            turnAllowed;
1841    const AlignParams& ali_params;
1842    int                maxProtection;         // protection level
1843
1844    // -------------------- new members
1845    int                wasNotAllowedToAlign;  // number of failures caused by wrong protection
1846    int                err_count;             // count errors
1847    bool               continue_on_error;         /* true -> run single alignments in separate transactions.
1848                                               *         If one target fails, continue with rest.
1849                                               * false -> run all in one transaction
1850                                               *          One fails -> all fail!
1851                                               */
1852    FA_errorAction     error_action;
1853
1854    typedef std::list<GBDATA*> GBDATAlist;
1855    GBDATAlist species_to_mark;       // species that will be marked after aligning
1856
1857    ARB_ERROR alignToReference(GBDATA *gb_toalign, const AlignmentReference& ref);
1858    ARB_ERROR alignTargetsToReference(const AlignmentReference& ref, GBDATA *gb_species_data);
1859
1860    ARB_ERROR alignToExplicitReference(GBDATA *gb_species_data, int max_seq_length);
1861    ARB_ERROR alignToConsensus(GBDATA *gb_species_data, int max_seq_length);
1862    ARB_ERROR alignToRelatives(GBDATA *gb_species_data, int max_seq_length);
1863
1864    void triggerAction(GBDATA *gb_species, bool has_been_aligned) {
1865        bool mark = false;
1866        switch (error_action) {
1867            case FA_MARK_FAILED:  mark = !has_been_aligned; break;
1868            case FA_MARK_ALIGNED: mark = has_been_aligned; break;
1869            case FA_NO_ACTION:    mark = false; break;
1870        }
1871        if (mark) species_to_mark.push_back(gb_species);
1872    }
1873
1874public:
1875
1876#if defined(WARN_TODO)
1877#warning pass AlignmentReference from caller (replacing reference parameters)
1878#endif
1879
1880    Aligner(GBDATA *gb_main_,
1881
1882            // define alignment target(s):
1883            FA_alignTarget                     alignWhat_,
1884            GB_CSTR                            alignment_,                   // name of alignment to use (if NULp -> use default)
1885            GB_CSTR                            toalign_,                     // name of species to align (used if alignWhat == FA_CURRENT)
1886            Aligner_get_first_selected_species get_first_selected_species_,  // used if alignWhat == FA_SELECTED
1887            Aligner_get_next_selected_species  get_next_selected_species_,   // --- "" ---
1888
1889            // define reference sequence(s):
1890            GB_CSTR                    reference_,     // name of reference species
1891            Aligner_get_consensus_func get_consensus_, // function to get consensus seq
1892            SearchRelativeParams&      relSearch_,     // params to search for relatives
1893
1894            // general params:
1895            FA_turn            turnAllowed_,
1896            const AlignParams& ali_params_,
1897            int                maxProtection_,      // protection level
1898            bool               continue_on_error_,
1899            FA_errorAction     error_action_)
1900        : gb_main(gb_main_),
1901          alignWhat(alignWhat_),
1902          alignment(alignment_),
1903          toalign(toalign_),
1904          get_first_selected_species(get_first_selected_species_),
1905          get_next_selected_species(get_next_selected_species_),
1906          reference(reference_),
1907          get_consensus(get_consensus_),
1908          relSearch(relSearch_),
1909          turnAllowed(turnAllowed_),
1910          ali_params(ali_params_),
1911          maxProtection(maxProtection_),
1912          wasNotAllowedToAlign(0),
1913          err_count(0),
1914          continue_on_error(continue_on_error_),
1915          error_action(continue_on_error ? error_action_ : FA_NO_ACTION)
1916    {}
1917
1918    ARB_ERROR run();
1919};
1920
1921ARB_ERROR Aligner::alignToReference(GBDATA *gb_toalign, const AlignmentReference& ref) {
1922    int       myProtection = GB_read_security_write(GBT_find_sequence(gb_toalign, alignment));
1923    ARB_ERROR error;
1924
1925    if (myProtection<=maxProtection) {
1926        // Depending on 'continue_on_error' we either
1927        // * stop aligning if an error occurs or
1928        // * run the alignment of each species in its own transaction, ignore but report errors and
1929        //   optionally mark aligned or failed species.
1930       
1931        if (continue_on_error) {
1932            fa_assert(GB_get_transaction_level(gb_main) == 1);
1933            error = GB_end_transaction(gb_main, error); // end global transaction
1934        }
1935
1936        if (!error) {
1937            error             = GB_push_transaction(gb_main);
1938            if (!error) error = ref.align_to(gb_toalign);
1939            error             = GB_end_transaction(gb_main, error);
1940
1941            if (error) err_count++;
1942            triggerAction(gb_toalign, !error);
1943        }
1944
1945        if (continue_on_error) {
1946            if (error) {
1947                GB_warning(error.deliver());
1948                error = NULp;
1949            }
1950            error = GB_begin_transaction(gb_main); // re-open global transaction
1951        }
1952    }
1953    else {
1954        wasNotAllowedToAlign++;
1955        triggerAction(gb_toalign, false);
1956    }
1957
1958    return error;
1959}
1960
1961ARB_ERROR Aligner::alignTargetsToReference(const AlignmentReference& ref, GBDATA *gb_species_data) {
1962    ARB_ERROR error;
1963
1964    fa_assert(GB_get_transaction_level(gb_main) == 1);
1965   
1966    switch (alignWhat) {
1967        case FA_CURRENT: { // align one sequence
1968            fa_assert(toalign);
1969
1970            GBDATA *gb_toalign = GBT_find_species_rel_species_data(gb_species_data, toalign);
1971            if (!gb_toalign) {
1972                error = species_not_found(toalign);
1973            }
1974            else {
1975                currentSequenceNumber = overallSequenceNumber = 1;
1976                error = alignToReference(gb_toalign, ref);
1977            }
1978            break;
1979        }
1980        case FA_MARKED: { // align all marked sequences
1981            int     count      = GBT_count_marked_species(gb_main);
1982            GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data);
1983
1984            arb_progress progress("Aligning marked species", long(count));
1985            progress.auto_subtitles("Species");
1986
1987            currentSequenceNumber = 1;
1988            overallSequenceNumber = count;
1989
1990            while (gb_species && !error) {
1991                error      = alignToReference(gb_species, ref);
1992                progress.inc_and_check_user_abort(error);
1993                gb_species = GBT_next_marked_species(gb_species);
1994            }
1995            break;
1996        }
1997        case FA_SELECTED: { // align all selected species
1998            int     count;
1999            GBDATA *gb_species = get_first_selected_species(&count);
2000
2001           
2002            currentSequenceNumber = 1;
2003            overallSequenceNumber = count;
2004
2005            if (!gb_species) {
2006                aw_message("There is no selected species!");
2007            }
2008            else {
2009                arb_progress progress("Aligning selected species", long(count));
2010                progress.auto_subtitles("Species");
2011
2012                while (gb_species && !error) {
2013                    error      = alignToReference(gb_species, ref);
2014                    progress.inc_and_check_user_abort(error);
2015                    gb_species = get_next_selected_species();
2016                }
2017            }
2018            break;
2019        }
2020    }
2021   
2022    fa_assert(GB_get_transaction_level(gb_main) == 1);
2023    return error;
2024}
2025
2026ARB_ERROR Aligner::alignToExplicitReference(GBDATA *gb_species_data, int max_seq_length) {
2027    ARB_ERROR  error;
2028    GBDATA    *gb_reference = GBT_find_species_rel_species_data(gb_species_data, reference);
2029
2030    if (!gb_reference) {
2031        error = species_not_found(reference);
2032    }
2033    else {
2034        long                  referenceChksum;
2035        CompactedSubSequence *referenceSeq = readCompactedSequence(gb_reference, alignment, &error, NULp, &referenceChksum, ali_params.range);
2036
2037        if (island_hopper) {
2038#if defined(WARN_TODO)
2039#warning setting island_hopper reference has to be done in called function (seems that it is NOT done for alignToConsensus and alignToRelatives). First get tests in place!
2040#endif
2041            GBDATA *gb_seq = GBT_find_sequence(gb_reference, alignment);        // get sequence
2042            if (gb_seq) {
2043                long        length = GB_read_string_count(gb_seq);
2044                const char *data   = GB_read_char_pntr(gb_seq);
2045
2046                island_hopper->set_ref_sequence(data);
2047                island_hopper->set_alignment_length(length);
2048            }
2049        }
2050
2051
2052        if (!error) {
2053#if defined(WARN_TODO)
2054#warning do not pass FastSearchSequence to ExplicitReference, instead pass sequence and length (ExplicitReference shall create it itself)
2055#endif
2056
2057            FastSearchSequence referenceFastSeq(*referenceSeq);
2058            ExplicitReference  target(alignment, &referenceFastSeq, gb_reference, max_seq_length, ali_params);
2059
2060            error = alignTargetsToReference(target, gb_species_data);
2061        }
2062        delete referenceSeq;
2063    }
2064    return error;
2065}
2066
2067ARB_ERROR Aligner::alignToConsensus(GBDATA *gb_species_data, int max_seq_length) {
2068    return alignTargetsToReference(ConsensusReference(alignment, get_consensus, max_seq_length, ali_params),
2069                                   gb_species_data);
2070}
2071
2072ARB_ERROR Aligner::alignToRelatives(GBDATA *gb_species_data, int max_seq_length) {
2073   
2074    return alignTargetsToReference(SearchRelativesReference(relSearch, max_seq_length, turnAllowed, alignment, ali_params),
2075                                   gb_species_data);
2076}
2077
2078ARB_ERROR Aligner::run() {
2079    fa_assert(GB_get_transaction_level(gb_main) == 0); // no open transaction allowed here!
2080    ARB_ERROR error = GB_push_transaction(gb_main);
2081
2082    // if neither 'reference' nor 'get_consensus' are specified -> use next relative for (each) sequence:
2083    bool search_by_pt_server = !reference && !get_consensus;
2084
2085    err_count            = 0;
2086    wasNotAllowedToAlign = 0;                       // incremented for every sequence which has higher protection level (and was not aligned)
2087    species_to_mark.clear();
2088
2089    fa_assert(!reference || !get_consensus);    // can't do both modes
2090
2091    if (turnAllowed != FA_TURN_NEVER) {
2092        if ((ali_params.range.is_part()) || !search_by_pt_server) {
2093            // if not selected 'Range/Whole sequence' or not selected 'Reference/Auto search..'
2094            turnAllowed = FA_TURN_NEVER; // then disable mirroring for the current call
2095        }
2096    }
2097
2098    if (!error && !alignment) {
2099        alignment = (GB_CSTR)GBT_get_default_alignment(gb_main); // get default alignment
2100        if (!alignment) error = "No default alignment";
2101    }
2102
2103    if (!error && alignment) {
2104        global_alignmentType = GBT_get_alignment_type(gb_main, alignment);
2105        if (search_by_pt_server) {
2106            GB_alignment_type pt_server_alignmentType = GBT_get_alignment_type(gb_main, relSearch.pt_server_alignment);
2107
2108            if (pt_server_alignmentType != GB_AT_RNA &&
2109                pt_server_alignmentType != GB_AT_DNA) {
2110                error = "pt_servers only support RNA/DNA sequences.\n"
2111                    "In the aligner window you may specify a RNA/DNA alignment \n"
2112                    "and use a pt_server build on that alignment.";
2113            }
2114        }
2115    }
2116
2117    if (!error) {
2118        GBDATA *gb_species_data = GBT_get_species_data(gb_main);
2119        int max_seq_length = GBT_get_alignment_len(gb_main, alignment);
2120
2121        if (reference) error          = alignToExplicitReference(gb_species_data, max_seq_length);
2122        else if (get_consensus) error = alignToConsensus(gb_species_data, max_seq_length);
2123        else error                    = alignToRelatives(gb_species_data, max_seq_length);
2124    }
2125
2126    ClustalV_exit();
2127    unaligned_bases.clear();
2128
2129    error = GB_end_transaction(gb_main, error); // close global transaction
2130
2131    if (wasNotAllowedToAlign>0) {
2132        const char *mess = GBS_global_string("%i species were not aligned (because of protection level)", wasNotAllowedToAlign);
2133        aw_message(mess);
2134    }
2135
2136    if (err_count) {
2137        aw_message_if(error);
2138        error = GBS_global_string("Aligner produced %i error%c", err_count, err_count==1 ? '\0' : 's');
2139    }
2140
2141    if (error_action != FA_NO_ACTION) {
2142        fa_assert(continue_on_error);
2143
2144        GB_transaction ta(gb_main);
2145        GBT_mark_all(gb_main, 0);
2146        for (GBDATAlist::iterator sp = species_to_mark.begin(); sp != species_to_mark.end(); ++sp) {
2147            GB_write_flag(*sp, 1);
2148        }
2149
2150        const char *whatsMarked = (error_action == FA_MARK_ALIGNED) ? "aligned" : "failed";
2151        size_t      markCount   = species_to_mark.size();
2152        if (markCount>0) {
2153            const char *msg = GBS_global_string("%zu %s species %s been marked",
2154                                                markCount,
2155                                                whatsMarked,
2156                                                (markCount == 1) ? "has" : "have");
2157            aw_message(msg);
2158        }
2159    }
2160
2161    return error;
2162}
2163
2164void FastAligner_start(AW_window *aw, const AlignDataAccess *data_access) {
2165    AW_root   *root          = aw->get_root();
2166    char      *reference     = NULp; // align against next relatives
2167    char      *toalign       = NULp; // align marked species
2168    ARB_ERROR  error         = NULp;
2169    int        get_consensus = 0;
2170    int        pt_server_id  = -1;
2171
2172    Aligner_get_first_selected_species get_first_selected_species = NULp;
2173    Aligner_get_next_selected_species  get_next_selected_species  = NULp;
2174
2175    fa_assert(!island_hopper);
2176    if (root->awar(FA_AWAR_USE_ISLAND_HOPPING)->read_int()) {
2177        island_hopper = new IslandHopping;
2178        if (root->awar(FA_AWAR_USE_SECONDARY)->read_int()) {
2179            char *helix_string = data_access->getHelixString();
2180            if (helix_string) {
2181                island_hopper->set_helix(helix_string);
2182                free(helix_string);
2183            }
2184            else {
2185                error = "Warning: No HELIX found. Can't use secondary structure";
2186            }
2187        }
2188
2189        if (!error) {
2190            island_hopper->set_parameters(root->awar(FA_AWAR_ESTIMATE_BASE_FREQ)->read_int(),
2191                                          root->awar(FA_AWAR_BASE_FREQ_T)->read_float(),
2192                                          root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2193                                          root->awar(FA_AWAR_BASE_FREQ_A)->read_float(),
2194                                          root->awar(FA_AWAR_BASE_FREQ_C)->read_float(),
2195                                          root->awar(FA_AWAR_SUBST_PARA_CT)->read_float(),
2196                                          root->awar(FA_AWAR_SUBST_PARA_AT)->read_float(),
2197                                          root->awar(FA_AWAR_SUBST_PARA_GT)->read_float(),
2198                                          root->awar(FA_AWAR_SUBST_PARA_AC)->read_float(),
2199                                          root->awar(FA_AWAR_SUBST_PARA_CG)->read_float(),
2200                                          root->awar(FA_AWAR_SUBST_PARA_AG)->read_float(),
2201                                          root->awar(FA_AWAR_EXPECTED_DISTANCE)->read_float(),
2202                                          root->awar(FA_AWAR_STRUCTURE_SUPPLEMENT)->read_float(),
2203                                          root->awar(FA_AWAR_GAP_A)->read_float(),
2204                                          root->awar(FA_AWAR_GAP_B)->read_float(),
2205                                          root->awar(FA_AWAR_GAP_C)->read_float(),
2206                                          root->awar(FA_AWAR_THRESHOLD)->read_float()
2207                                          );
2208        }
2209    }
2210
2211    FA_alignTarget alignWhat = static_cast<FA_alignTarget>(root->awar(FA_AWAR_TO_ALIGN)->read_int());
2212    if (!error) {
2213        switch (alignWhat) {
2214            case FA_CURRENT: { // align current species
2215                toalign = root->awar(AWAR_SPECIES_NAME)->read_string();
2216                break;
2217            }
2218            case FA_MARKED: { // align marked species
2219                break;
2220            }
2221            case FA_SELECTED: { // align selected species
2222                get_first_selected_species = data_access->get_first_selected_species;
2223                get_next_selected_species  = data_access->get_next_selected_species;
2224                break;
2225            }
2226            default: {
2227                fa_assert(0);
2228                break;
2229            }
2230        }
2231
2232        switch (static_cast<FA_reference>(root->awar(FA_AWAR_REFERENCE)->read_int())) {
2233            case FA_REF_EXPLICIT: // align against specified species
2234                reference = root->awar(FA_AWAR_REFERENCE_NAME)->read_string();
2235                break;
2236
2237            case FA_REF_CONSENSUS: // align against group consensus
2238                if (data_access->get_group_consensus) {
2239                    get_consensus = 1;
2240                }
2241                else {
2242                    error = "Can't get group consensus here.";
2243                }
2244                break;
2245
2246            case FA_REF_RELATIVES: // align against species searched via pt_server
2247                pt_server_id = root->awar(AWAR_PT_SERVER)->read_int();
2248                if (pt_server_id<0) {
2249                    error = "No pt_server selected";
2250                }
2251                break;
2252
2253            default: fa_assert(0);
2254                break;
2255        }
2256    }
2257
2258    RangeList ranges;
2259    bool      autoRestrictRange4nextRelSearch = true;
2260
2261    if (!error) {
2262        switch (static_cast<FA_range>(root->awar(FA_AWAR_RANGE)->read_int())) {
2263            case FA_WHOLE_SEQUENCE:
2264                autoRestrictRange4nextRelSearch = false;
2265                ranges.add(PosRange::whole());
2266                break;
2267
2268            case FA_AROUND_CURSOR: {
2269                int curpos = root->awar(AWAR_CURSOR_POSITION_LOCAL)->read_int();
2270                int size = root->awar(FA_AWAR_AROUND)->read_int();
2271
2272                ranges.add(PosRange(curpos-size, curpos+size));
2273                break;
2274            }
2275            case FA_SELECTED_RANGE: {
2276                PosRange range;
2277                if (!data_access->get_selected_range(range)) {
2278                    error = "There is no selected species!";
2279                }
2280                else {
2281                    ranges.add(range);
2282                }
2283                break;
2284            }
2285
2286            case FA_SAI_MULTI_RANGE: {
2287                GB_transaction ta(data_access->gb_main);
2288
2289                char *sai_name = root->awar(FA_AWAR_SAI_RANGE_NAME)->read_string();
2290                char *aliuse   = GBT_get_default_alignment(data_access->gb_main);
2291
2292                GBDATA *gb_sai     = GBT_expect_SAI(data_access->gb_main, sai_name);
2293                if (!gb_sai) error = GB_await_error();
2294                else {
2295                    GBDATA *gb_data = GBT_find_sequence(gb_sai, aliuse);
2296                    if (!gb_data) {
2297                        error = GB_have_error()
2298                            ? GB_await_error()
2299                            : GBS_global_string("SAI '%s' has no data in alignment '%s'", sai_name, aliuse);
2300                    }
2301                    else {
2302                        char *sai_data = GB_read_string(gb_data);
2303                        char *set_bits = root->awar(FA_AWAR_SAI_RANGE_CHARS)->read_string();
2304
2305                        ranges = build_RangeList_from_string(sai_data, set_bits, false);
2306
2307                        free(set_bits);
2308                        free(sai_data);
2309                    }
2310                }
2311                free(aliuse);
2312                free(sai_name);
2313                break;
2314            }
2315        }
2316    }
2317
2318    if (!error) {
2319        for (RangeList::iterator r = ranges.begin(); r != ranges.end() && !error; ++r) {
2320            PosRange range = *r;
2321
2322            GBDATA *gb_main          = data_access->gb_main;
2323            char   *editor_alignment = NULp;
2324            long    alignment_length;
2325            {
2326                GB_transaction  ta(gb_main);
2327                char           *default_alignment = GBT_get_default_alignment(gb_main);
2328
2329                alignment_length = GBT_get_alignment_len(gb_main, default_alignment);
2330
2331                editor_alignment = root->awar_string(AWAR_EDITOR_ALIGNMENT, default_alignment)->read_string();
2332                free(default_alignment);
2333
2334            }
2335
2336            char     *pt_server_alignment = root->awar(FA_AWAR_PT_SERVER_ALIGNMENT)->read_string();
2337            PosRange  relRange            = PosRange::whole(); // unrestricted!
2338
2339            if (autoRestrictRange4nextRelSearch) {
2340                AW_awar    *awar_relrange = root->awar(FA_AWAR_RELATIVE_RANGE);
2341                const char *relrange      = awar_relrange->read_char_pntr();
2342                if (relrange[0]) {
2343                    int region_plus = atoi(relrange);
2344
2345                    fa_assert(range.is_part());
2346
2347                    relRange = PosRange(range.start()-region_plus, range.end()+region_plus); // restricted
2348                    awar_relrange->write_as_string(GBS_global_string("%i", region_plus)); // set awar to detected value (avoid misbehavior when it contains ' ' or similar)
2349                }
2350            }
2351
2352            if (island_hopper) {
2353                island_hopper->set_range(ExplicitRange(range, alignment_length));
2354                range = PosRange::whole();
2355            }
2356
2357            SearchRelativeParams relSearch(new PT_FamilyFinder(gb_main,
2358                                                               pt_server_id,
2359                                                               root->awar(AWAR_NN_OLIGO_LEN)->read_int(),
2360                                                               root->awar(AWAR_NN_MISMATCHES)->read_int(),
2361                                                               root->awar(AWAR_NN_FAST_MODE)->read_int(),
2362                                                               root->awar(AWAR_NN_REL_MATCHES)->read_int(),
2363                                                               RSS_BOTH_MIN), // old scaling as b4 [8520] @@@ make configurable
2364                                           pt_server_alignment,
2365                                           root->awar(FA_AWAR_NEXT_RELATIVES)->read_int());
2366
2367            relSearch.getFamilyFinder()->restrict_2_region(relRange);
2368
2369            struct AlignParams ali_params = {
2370                static_cast<FA_report>(root->awar(FA_AWAR_REPORT)->read_int()),
2371                bool(root->awar(FA_AWAR_SHOW_GAPS_MESSAGES)->read_int()),
2372                range
2373            };
2374
2375            {
2376                arb_progress progress("FastAligner");
2377                progress.allow_title_reuse();
2378
2379                int cont_on_error = root->awar(FA_AWAR_CONTINUE_ON_ERROR)->read_int();
2380
2381                Aligner aligner(gb_main,
2382                                alignWhat,
2383                                editor_alignment,
2384                                toalign,
2385                                get_first_selected_species,
2386                                get_next_selected_species,
2387
2388                                reference,
2389                                get_consensus ? data_access->get_group_consensus : NULp,
2390                                relSearch,
2391
2392                                static_cast<FA_turn>(root->awar(FA_AWAR_MIRROR)->read_int()),
2393                                ali_params,
2394                                root->awar(FA_AWAR_PROTECTION)->read_int(),
2395                                cont_on_error,
2396                                (FA_errorAction)root->awar(FA_AWAR_ACTION_ON_ERROR)->read_int());
2397                error = aligner.run();
2398
2399                if (error && cont_on_error) {
2400                    aw_message_if(error);
2401                    error = NULp;
2402                }
2403            }
2404
2405            free(pt_server_alignment);
2406            free(editor_alignment);
2407        }
2408    }
2409
2410    if (island_hopper) {
2411        delete island_hopper;
2412        island_hopper = NULp;
2413    }
2414
2415    if (toalign) free(toalign);
2416    aw_message_if(error);
2417    if (data_access->do_refresh) data_access->refresh_display();
2418}
2419
2420static void nullcb() { }
2421
2422void FastAligner_create_variables(AW_root *root, AW_default db1) {
2423    root->awar_string(FA_AWAR_REFERENCE_NAME, "", db1);
2424
2425    root->awar_int(FA_AWAR_TO_ALIGN,  FA_CURRENT,        db1);
2426    root->awar_int(FA_AWAR_REFERENCE, FA_REF_EXPLICIT,   db1);
2427    root->awar_int(FA_AWAR_RANGE,     FA_WHOLE_SEQUENCE, db1);
2428
2429    AW_awar *ali_protect = root->awar_int(FA_AWAR_PROTECTION, 0, db1);
2430    if (ARB_in_novice_mode(root)) {
2431        ali_protect->write_int(0); // reset protection for noobs
2432    }
2433
2434    root->awar_int(FA_AWAR_AROUND,             25,                  db1);
2435    root->awar_int(FA_AWAR_MIRROR,             FA_TURN_INTERACTIVE, db1);
2436    root->awar_int(FA_AWAR_REPORT,             FA_NO_REPORT,        db1);
2437    root->awar_int(FA_AWAR_SHOW_GAPS_MESSAGES, 1,                   db1);
2438    root->awar_int(FA_AWAR_CONTINUE_ON_ERROR,  1,                   db1);
2439    root->awar_int(FA_AWAR_ACTION_ON_ERROR,    FA_NO_ACTION,        db1);
2440    root->awar_int(FA_AWAR_USE_SECONDARY,      0,                   db1);
2441    root->awar_int(AWAR_PT_SERVER,             0,                   db1);
2442    root->awar_int(FA_AWAR_NEXT_RELATIVES,     1,                   db1)->set_minmax(1, 100);
2443
2444    root->awar_string(FA_AWAR_RELATIVE_RANGE, "", db1);
2445    root->awar_string(FA_AWAR_PT_SERVER_ALIGNMENT, root->awar(AWAR_DEFAULT_ALIGNMENT)->read_char_pntr(), db1);
2446
2447    root->awar_string(FA_AWAR_SAI_RANGE_NAME,  "",   db1);
2448    root->awar_string(FA_AWAR_SAI_RANGE_CHARS, "x1", db1);
2449
2450    // island hopping:
2451
2452    root->awar_int(FA_AWAR_USE_ISLAND_HOPPING, 0, db1);
2453
2454    root->awar_int(FA_AWAR_ESTIMATE_BASE_FREQ, 1, db1);
2455
2456    root->awar_float(FA_AWAR_BASE_FREQ_A, 0.25, db1);
2457    root->awar_float(FA_AWAR_BASE_FREQ_C, 0.25, db1);
2458    root->awar_float(FA_AWAR_BASE_FREQ_G, 0.25, db1);
2459    root->awar_float(FA_AWAR_BASE_FREQ_T, 0.25, db1);
2460
2461    root->awar_float(FA_AWAR_SUBST_PARA_AC, 1.0, db1);
2462    root->awar_float(FA_AWAR_SUBST_PARA_AG, 4.0, db1);
2463    root->awar_float(FA_AWAR_SUBST_PARA_AT, 1.0, db1);
2464    root->awar_float(FA_AWAR_SUBST_PARA_CG, 1.0, db1);
2465    root->awar_float(FA_AWAR_SUBST_PARA_CT, 4.0, db1);
2466    root->awar_float(FA_AWAR_SUBST_PARA_GT, 1.0, db1);
2467
2468    root->awar_float(FA_AWAR_EXPECTED_DISTANCE,    0.3,   db1);
2469    root->awar_float(FA_AWAR_STRUCTURE_SUPPLEMENT, 0.5,   db1);
2470    root->awar_float(FA_AWAR_THRESHOLD,            0.005, db1);
2471
2472    root->awar_float(FA_AWAR_GAP_A, 8.0, db1);
2473    root->awar_float(FA_AWAR_GAP_B, 4.0, db1);
2474    root->awar_float(FA_AWAR_GAP_C, 7.0, db1);
2475
2476    AWTC_create_common_next_neighbour_vars(root, makeRootCallback(nullcb));
2477}
2478
2479void FastAligner_set_align_current(AW_root *root, AW_default db1) {
2480    root->awar_int(FA_AWAR_TO_ALIGN, FA_CURRENT, db1);
2481}
2482
2483void FastAligner_set_reference_species(AW_root *root) {
2484    // sets the aligner reference species to current species
2485    char *specName = root->awar(AWAR_SPECIES_NAME)->read_string();
2486    root->awar(FA_AWAR_REFERENCE_NAME)->write_string(specName);
2487    free(specName);
2488}
2489
2490static AW_window *create_island_hopping_window(AW_root *root) {
2491    AW_window_simple *aws = new AW_window_simple;
2492
2493    aws->init(root, "ISLAND_HOPPING_PARA", "Parameters for Island Hopping");
2494    aws->load_xfig("faligner/islandhopping.fig");
2495
2496    aws->at("close");
2497    aws->callback(AW_POPDOWN);
2498    aws->create_button("CLOSE", "CLOSE", "O");
2499
2500    aws->at("help");
2501    aws->callback(makeHelpCallback("islandhopping.hlp"));
2502    aws->create_button("HELP", "HELP");
2503
2504    aws->at("use_secondary");
2505    aws->label("Use secondary structure (only for re-align)");
2506    aws->create_toggle(FA_AWAR_USE_SECONDARY);
2507
2508    aws->at("freq");
2509    aws->create_toggle_field(FA_AWAR_ESTIMATE_BASE_FREQ, "Base freq.", "B");
2510    aws->insert_default_toggle("Estimate", "E", 1);
2511    aws->insert_toggle("Define here: ", "D", 0);
2512    aws->update_toggle_field();
2513
2514    aws->at("freq_a"); aws->label("A:"); aws->create_input_field(FA_AWAR_BASE_FREQ_A, 4);
2515    aws->at("freq_c"); aws->label("C:"); aws->create_input_field(FA_AWAR_BASE_FREQ_C, 4);
2516    aws->at("freq_g"); aws->label("G:"); aws->create_input_field(FA_AWAR_BASE_FREQ_G, 4);
2517    aws->at("freq_t"); aws->label("T:"); aws->create_input_field(FA_AWAR_BASE_FREQ_T, 4);
2518
2519    int xpos[4], ypos[4];
2520    {
2521        aws->button_length(1);
2522
2523        int dummy;
2524        aws->at("h_a"); aws->get_at_position(&xpos[0], &dummy); aws->create_button(NULp, "A");
2525        aws->at("h_c"); aws->get_at_position(&xpos[1], &dummy); aws->create_button(NULp, "C");
2526        aws->at("h_g"); aws->get_at_position(&xpos[2], &dummy); aws->create_button(NULp, "G");
2527        aws->at("h_t"); aws->get_at_position(&xpos[3], &dummy); aws->create_button(NULp, "T");
2528
2529        aws->at("v_a"); aws->get_at_position(&dummy, &ypos[0]); aws->create_button(NULp, "A");
2530        aws->at("v_c"); aws->get_at_position(&dummy, &ypos[1]); aws->create_button(NULp, "C");
2531        aws->at("v_g"); aws->get_at_position(&dummy, &ypos[2]); aws->create_button(NULp, "G");
2532        aws->at("v_t"); aws->get_at_position(&dummy, &ypos[3]); aws->create_button(NULp, "T");
2533    }
2534
2535    aws->at("subst"); aws->create_button(NULp, "Substitution rate parameters:");
2536
2537#define XOFF -25
2538#define YOFF 0
2539
2540    aws->at(xpos[1]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AC, 4);
2541    aws->at(xpos[2]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AG, 4);
2542    aws->at(xpos[3]+XOFF, ypos[0]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_AT, 4);
2543    aws->at(xpos[2]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CG, 4);
2544    aws->at(xpos[3]+XOFF, ypos[1]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_CT, 4);
2545    aws->at(xpos[3]+XOFF, ypos[2]+YOFF); aws->create_input_field(FA_AWAR_SUBST_PARA_GT, 4);
2546
2547#undef XOFF
2548#undef YOFF
2549
2550    aws->label_length(22);
2551
2552    aws->at("dist");
2553    aws->label("Expected distance");
2554    aws->create_input_field(FA_AWAR_EXPECTED_DISTANCE, 5);
2555
2556    aws->at("supp");
2557    aws->label("Structure supplement");
2558    aws->create_input_field(FA_AWAR_STRUCTURE_SUPPLEMENT, 5);
2559
2560    aws->at("thres");
2561    aws->label("Threshold");
2562    aws->create_input_field(FA_AWAR_THRESHOLD, 5);
2563
2564    aws->label_length(10);
2565
2566    aws->at("gapA");
2567    aws->label("Gap A");
2568    aws->create_input_field(FA_AWAR_GAP_A, 5);
2569
2570    aws->at("gapB");
2571    aws->label("Gap B");
2572    aws->create_input_field(FA_AWAR_GAP_B, 5);
2573
2574    aws->at("gapC");
2575    aws->label("Gap C");
2576    aws->create_input_field(FA_AWAR_GAP_C, 5);
2577
2578    return aws;
2579}
2580
2581static AW_window *create_family_settings_window(AW_root *root) {
2582    static AW_window_simple *aws = NULp;
2583
2584    if (!aws) {
2585        aws = new AW_window_simple;
2586
2587        aws->init(root, "FAMILY_PARAMS", "Family search parameters");
2588        aws->load_xfig("faligner/family_settings.fig");
2589
2590        aws->at("close");
2591        aws->callback(AW_POPDOWN);
2592        aws->create_button("CLOSE", "CLOSE", "O");
2593
2594        aws->at("help");
2595        aws->callback(makeHelpCallback("next_neighbours_common.hlp"));
2596        aws->create_button("HELP", "HELP");
2597
2598        aws->auto_space(5, 5);
2599        AWTC_create_common_next_neighbour_fields(aws, 200);
2600    }
2601
2602    return aws;
2603}
2604
2605static AWT_config_mapping_def aligner_config_mapping[] = {
2606    { FA_AWAR_USE_ISLAND_HOPPING,  "island" },
2607    { FA_AWAR_TO_ALIGN,            "target" },
2608    { FA_AWAR_REFERENCE,           "ref" },
2609    { FA_AWAR_REFERENCE_NAME,      "refname" },
2610    { FA_AWAR_RELATIVE_RANGE,      "relrange" },
2611    { FA_AWAR_NEXT_RELATIVES,      "relatives" },
2612    { FA_AWAR_PT_SERVER_ALIGNMENT, "ptali" },
2613
2614    // relative-search specific parameters from subwindow (create_family_settings_window)
2615    // same as ../DB_UI/ui_species.cxx@RELATIVES_CONFIG
2616    { AWAR_NN_OLIGO_LEN,   "oligolen" },
2617    { AWAR_NN_MISMATCHES,  "mismatches" },
2618    { AWAR_NN_FAST_MODE,   "fastmode" },
2619    { AWAR_NN_REL_MATCHES, "relmatches" },
2620    { AWAR_NN_REL_SCALING, "relscaling" },
2621
2622    // island-hopping parameters (create_island_hopping_window)
2623    { FA_AWAR_USE_SECONDARY,        "use2nd" },
2624    { FA_AWAR_ESTIMATE_BASE_FREQ,   "estbasefreq" },
2625    { FA_AWAR_BASE_FREQ_A,          "freq_a" },
2626    { FA_AWAR_BASE_FREQ_C,          "freq_c" },
2627    { FA_AWAR_BASE_FREQ_G,          "freq_g" },
2628    { FA_AWAR_BASE_FREQ_T,          "freq_t" },
2629    { FA_AWAR_SUBST_PARA_AC,        "subst_ac" },
2630    { FA_AWAR_SUBST_PARA_AG,        "subst_ag" },
2631    { FA_AWAR_SUBST_PARA_AT,        "subst_at" },
2632    { FA_AWAR_SUBST_PARA_CG,        "subst_cg" },
2633    { FA_AWAR_SUBST_PARA_CT,        "subst_ct" },
2634    { FA_AWAR_SUBST_PARA_GT,        "subst_gt" },
2635    { FA_AWAR_EXPECTED_DISTANCE,    "distance" },
2636    { FA_AWAR_STRUCTURE_SUPPLEMENT, "supplement" },
2637    { FA_AWAR_THRESHOLD,            "threshold" },
2638    { FA_AWAR_GAP_A,                "gap_a" },
2639    { FA_AWAR_GAP_B,                "gap_b" },
2640    { FA_AWAR_GAP_C,                "gap_c" },
2641
2642    { NULp, NULp }
2643};
2644
2645AW_window *FastAligner_create_window(AW_root *root, const AlignDataAccess *data_access) {
2646    AW_window_simple *aws = new AW_window_simple;
2647
2648    aws->init(root, "INTEGRATED_ALIGNERS", INTEGRATED_ALIGNERS_TITLE);
2649    aws->load_xfig("faligner/faligner.fig");
2650
2651    aws->label_length(10);
2652    aws->button_length(10);
2653
2654    aws->at("close");
2655    aws->callback(AW_POPDOWN);
2656    aws->create_button("CLOSE", "CLOSE", "O");
2657
2658    aws->at("help");
2659    aws->callback(makeHelpCallback("faligner.hlp"));
2660    aws->create_button("HELP", "HELP");
2661
2662    aws->at("aligner");
2663    aws->create_toggle_field(FA_AWAR_USE_ISLAND_HOPPING, "Aligner", "A");
2664    aws->insert_default_toggle("Fast aligner",   "F", 0);
2665    aws->sens_mask(AWM_EXP);
2666    aws->insert_toggle        ("Island Hopping", "I", 1);
2667    aws->sens_mask(AWM_ALL);
2668    aws->update_toggle_field();
2669
2670    aws->button_length(12);
2671    aws->at("island_para");
2672    aws->callback(create_island_hopping_window);
2673    aws->sens_mask(AWM_EXP);
2674    aws->create_button("island_para", "Parameters", "");
2675    aws->sens_mask(AWM_ALL);
2676
2677    aws->button_length(10);
2678
2679    aws->at("rev_compl");
2680    aws->callback(makeWindowCallback(build_reverse_complement, data_access));
2681    aws->create_button("reverse_complement", "Turn now!", "");
2682
2683    aws->at("what");
2684    aws->create_toggle_field(FA_AWAR_TO_ALIGN, "Align what?", "A");
2685    aws->insert_toggle        ("Current Species:", "A", FA_CURRENT);
2686    aws->insert_default_toggle("Marked Species",   "M", FA_MARKED);
2687    aws->insert_toggle        ("Selected Species", "S", FA_SELECTED);
2688    aws->update_toggle_field();
2689
2690    aws->at("swhat");
2691    aws->create_input_field(AWAR_SPECIES_NAME, 2);
2692
2693    aws->at("against");
2694    aws->create_toggle_field(FA_AWAR_REFERENCE, "Reference", "");
2695    aws->insert_toggle        ("Species by name:",          "S", FA_REF_EXPLICIT);
2696    aws->insert_toggle        ("Group consensus",           "K", FA_REF_CONSENSUS);
2697    aws->insert_default_toggle("Auto search by pt_server:", "A", FA_REF_RELATIVES);
2698    aws->update_toggle_field();
2699
2700    aws->at("sagainst");
2701    aws->create_input_field(FA_AWAR_REFERENCE_NAME, 2);
2702
2703    aws->at("copy");
2704    aws->callback(RootAsWindowCallback::simple(FastAligner_set_reference_species));
2705    aws->create_button("Copy", "Copy", "");
2706
2707    aws->label_length(0);
2708    aws->at("pt_server");
2709    awt_create_PTSERVER_selection_button(aws, AWAR_PT_SERVER);
2710
2711    aws->label_length(23);
2712    aws->at("relrange");
2713    aws->label("Data from range only, plus");
2714    aws->create_input_field(FA_AWAR_RELATIVE_RANGE, 3);
2715   
2716    aws->at("relatives");
2717    aws->label("Number of relatives to use");
2718    aws->create_input_field(FA_AWAR_NEXT_RELATIVES, 3);
2719
2720    aws->label_length(9);
2721    aws->at("use_ali");
2722    aws->label("Alignment");
2723    aws->create_input_field(FA_AWAR_PT_SERVER_ALIGNMENT, 12);
2724
2725    aws->at("relSett");
2726    aws->callback(create_family_settings_window);
2727    aws->create_autosize_button("Settings", "More settings", "");
2728
2729    // Range
2730
2731    aws->label_length(10);
2732    aws->at("range");
2733    aws->create_toggle_field(FA_AWAR_RANGE, "Range", "");
2734    aws->insert_default_toggle("Whole sequence",            "", FA_WHOLE_SEQUENCE);
2735    aws->insert_toggle        ("Positions around cursor: ", "", FA_AROUND_CURSOR);
2736    aws->insert_toggle        ("Selected range",            "", FA_SELECTED_RANGE);
2737    aws->insert_toggle        ("Multi-Range by SAI",        "", FA_SAI_MULTI_RANGE);
2738    aws->update_toggle_field();
2739
2740    aws->at("around");
2741    aws->create_input_field(FA_AWAR_AROUND, 2);
2742
2743    aws->at("sai");
2744    awt_create_SAI_selection_button(data_access->gb_main, aws, FA_AWAR_SAI_RANGE_NAME);
2745   
2746    aws->at("rchars");
2747    aws->create_input_field(FA_AWAR_SAI_RANGE_CHARS, 2);
2748
2749    // Protection
2750
2751    aws->at("protection");
2752    aws->label("Protection");
2753    aws->create_option_menu(FA_AWAR_PROTECTION, true);
2754    aws->insert_default_option("0", NULp, 0);
2755    aws->insert_option        ("1", NULp, 1);
2756    aws->insert_option        ("2", NULp, 2);
2757    aws->insert_option        ("3", NULp, 3);
2758    aws->insert_option        ("4", NULp, 4);
2759    aws->insert_option        ("5", NULp, 5);
2760    aws->insert_option        ("6", NULp, 6);
2761    aws->update_option_menu();
2762
2763    // MirrorCheck
2764
2765    aws->at("mirror");
2766    aws->label("Turn check");
2767    aws->create_option_menu(FA_AWAR_MIRROR, true);
2768    aws->insert_option        ("Never turn sequence",         "", FA_TURN_NEVER);
2769    aws->insert_default_option("User acknowledgment ",        "", FA_TURN_INTERACTIVE);
2770    aws->insert_option        ("Automatically turn sequence", "", FA_TURN_ALWAYS);
2771    aws->update_option_menu();
2772
2773    // Report
2774
2775    aws->at("insert");
2776    aws->label("Report");
2777    aws->create_option_menu(FA_AWAR_REPORT, true);
2778    aws->insert_option        ("No report",                   "", FA_NO_REPORT);
2779    aws->sens_mask(AWM_EXP);
2780    aws->insert_default_option("Report to temporary entries", "", FA_TEMP_REPORT);
2781    aws->insert_option        ("Report to resident entries",  "", FA_REPORT);
2782    aws->sens_mask(AWM_ALL);
2783    aws->update_option_menu();
2784
2785    aws->at("gaps");
2786    aws->create_toggle(FA_AWAR_SHOW_GAPS_MESSAGES);
2787
2788    aws->at("continue");
2789    aws->create_toggle(FA_AWAR_CONTINUE_ON_ERROR);
2790
2791    aws->at("on_failure");
2792    aws->label("On failure");
2793    aws->create_option_menu(FA_AWAR_ACTION_ON_ERROR, true);
2794    aws->insert_default_option("do nothing",   "", FA_NO_ACTION);
2795    aws->insert_option        ("mark failed",  "", FA_MARK_FAILED);
2796    aws->insert_option        ("mark aligned", "", FA_MARK_ALIGNED);
2797    aws->update_option_menu();
2798
2799    aws->at("align");
2800    aws->callback(makeWindowCallback(FastAligner_start, data_access));
2801    aws->highlight();
2802    aws->create_button("GO", "GO", "G");
2803
2804    aws->at("config");
2805    AWT_insert_config_manager(aws, AW_ROOT_DEFAULT, "aligner", aligner_config_mapping);
2806
2807    return aws;
2808}
2809
2810// --------------------------------------------------------------------------------
2811
2812#ifdef UNIT_TESTS
2813
2814#include <test_unit.h>
2815
2816// ---------------------
2817//      OligoCounter
2818
2819#include <map>
2820#include <string>
2821
2822using std::map;
2823using std::string;
2824
2825typedef map<string, size_t> OligoCount;
2826
2827class OligoCounter {
2828    size_t oligo_len;
2829    size_t datasize;
2830   
2831    mutable OligoCount occurrence;
2832
2833    static string removeGaps(const char *seq) {
2834        size_t len = strlen(seq);
2835        string nogaps;
2836        nogaps.reserve(len);
2837
2838        for (size_t p = 0; p<len; ++p) {
2839            char c = seq[p];
2840            if (!is_gap(c)) nogaps.append(1, c);
2841        }
2842        return nogaps;
2843    }
2844
2845    void count_oligos(const string& seq) {
2846        occurrence.clear();
2847        size_t max_pos = seq.length()-oligo_len;
2848        for (size_t p = 0; p <= max_pos; ++p) {
2849            string oligo(seq, p, oligo_len);
2850            occurrence[oligo]++;
2851        }
2852    }
2853
2854public:
2855    OligoCounter()
2856        : oligo_len(0), 
2857          datasize(0)
2858    {}
2859    OligoCounter(const char *seq, size_t oligo_len_)
2860        : oligo_len(oligo_len_)
2861    {
2862        string seq_nogaps = removeGaps(seq);
2863        datasize          = seq_nogaps.length();
2864        count_oligos(seq_nogaps);
2865    }
2866
2867    size_t oligo_count(const char *oligo) {
2868        fa_assert(strlen(oligo) == oligo_len);
2869        return occurrence[oligo];
2870    }
2871
2872    size_t similarity_score(const OligoCounter& other) const {
2873        size_t score = 0;
2874        if (oligo_len == other.oligo_len) {
2875            for (OligoCount::const_iterator o = occurrence.begin(); o != occurrence.end(); ++o) {
2876                const string& oligo = o->first;
2877                size_t        count = o->second;
2878
2879                score += min(count, other.occurrence[oligo]);
2880            }
2881        }
2882        return score;
2883    }
2884
2885    size_t getDataSize() const { return datasize; }
2886};
2887
2888void TEST_OligoCounter() {
2889    OligoCounter oc1("CCAGGT", 3);
2890    OligoCounter oc2("GGTCCA", 3);
2891    OligoCounter oc2_gaps("..GGT--CCA..", 3);
2892    OligoCounter oc3("AGGTCC", 3);
2893    OligoCounter oc4("AGGTCCAGG", 3);
2894
2895    TEST_EXPECT_EQUAL(oc1.oligo_count("CCA"), 1);
2896    TEST_EXPECT_ZERO(oc1.oligo_count("CCG"));
2897    TEST_EXPECT_EQUAL(oc4.oligo_count("AGG"), 2);
2898
2899    int sc1_2 = oc1.similarity_score(oc2);
2900    int sc2_1 = oc2.similarity_score(oc1);
2901    TEST_EXPECT_EQUAL(sc1_2, sc2_1);
2902
2903    int sc1_2gaps = oc1.similarity_score(oc2_gaps);
2904    TEST_EXPECT_EQUAL(sc1_2, sc1_2gaps);
2905   
2906    int sc1_3 = oc1.similarity_score(oc3);
2907    int sc2_3 = oc2.similarity_score(oc3);
2908    int sc3_4 = oc3.similarity_score(oc4);
2909
2910    TEST_EXPECT_EQUAL(sc1_2, 2); // common oligos (CCA GGT)
2911    TEST_EXPECT_EQUAL(sc1_3, 2); // common oligos (AGG GGT)
2912    TEST_EXPECT_EQUAL(sc2_3, 3); // common oligos (GGT GTC TCC)
2913
2914    TEST_EXPECT_EQUAL(sc3_4, 4);
2915}
2916
2917// -------------------------
2918//      FakeFamilyFinder
2919
2920class FakeFamilyFinder: public FamilyFinder { // derived from a Noncopyable
2921    // used by unit tests to detect next relatives instead of asking the pt-server
2922
2923    GBDATA                    *gb_main;
2924    string                     ali_name;
2925    map<string, OligoCounter>  oligos_counted;      // key = species name
2926    PosRange                   counted_for_range;
2927    size_t                     oligo_len;
2928
2929public:
2930    FakeFamilyFinder(GBDATA *gb_main_, string ali_name_, bool rel_matches_, size_t oligo_len_)
2931        : FamilyFinder(rel_matches_, RSS_BOTH_MIN),
2932          gb_main(gb_main_),
2933          ali_name(ali_name_),
2934          counted_for_range(PosRange::whole()), 
2935          oligo_len(oligo_len_)
2936    {}
2937
2938    GB_ERROR searchFamily(const char *sequence, FF_complement compl_mode, int max_results, double min_score) OVERRIDE {
2939        // 'sequence' has to contain full sequence or part corresponding to 'range'
2940
2941        TEST_EXPECT_EQUAL(compl_mode, FF_FORWARD); // not fit for other modes
2942
2943        delete_family_list();
2944       
2945        OligoCounter seq_oligo_count(sequence, oligo_len);
2946
2947        if (range != counted_for_range) {
2948            oligos_counted.clear(); // forget results for different range
2949            counted_for_range = range;
2950        }
2951
2952        char *buffer     = NULp;
2953        int   buffersize = 0;
2954
2955        bool partial_match = range.is_part();
2956
2957        GB_transaction ta(gb_main);
2958        int            results = 0;
2959
2960        for (GBDATA *gb_species = GBT_first_species(gb_main);
2961             gb_species && results<max_results;
2962             gb_species = GBT_next_species(gb_species))
2963        {
2964            string name = GBT_get_name(gb_species);
2965            if (oligos_counted.find(name) == oligos_counted.end()) {
2966                GBDATA     *gb_data  = GBT_find_sequence(gb_species, ali_name.c_str());
2967                const char *spec_seq = GB_read_char_pntr(gb_data);
2968
2969                if (partial_match) {
2970                    int spec_seq_len = GB_read_count(gb_data);
2971                    int range_len    = ExplicitRange(range, spec_seq_len).size();
2972
2973                    if (buffersize<range_len) {
2974                        delete [] buffer;
2975                        buffersize = range_len;
2976                        buffer     = new char[buffersize+1];
2977                    }
2978
2979                    range.copy_corresponding_part(buffer, spec_seq, spec_seq_len);
2980                    oligos_counted[name] = OligoCounter(buffer, oligo_len);
2981                }
2982                else {
2983                    oligos_counted[name] = OligoCounter(spec_seq, oligo_len);
2984                }
2985            }
2986
2987            const OligoCounter& spec_oligo_count = oligos_counted[name];
2988            size_t              score            = seq_oligo_count.similarity_score(spec_oligo_count);
2989
2990            if (score>=min_score) {
2991                FamilyList *newMember = new FamilyList;
2992
2993                newMember->name        = strdup(name.c_str());
2994                newMember->matches     = score;
2995                newMember->rel_matches = score/spec_oligo_count.getDataSize();
2996                newMember->next        = NULp;
2997
2998                family_list = newMember->insertSortedBy_matches(family_list);
2999                results++;
3000            }
3001        }
3002
3003        delete [] buffer;
3004
3005        return NULp;
3006    }
3007};
3008
3009// ----------------------------
3010//      test_alignment_data
3011
3012#include <arb_unit_test.h>
3013
3014static const char *test_aliname = "ali_test";
3015
3016static const char *get_aligned_data_of(GBDATA *gb_main, const char *species_name) {
3017    GB_transaction  ta(gb_main);
3018    ARB_ERROR       error;
3019    const char     *data = NULp;
3020
3021    GBDATA *gb_species     = GBT_find_species(gb_main, species_name);
3022    if (!gb_species) error = GB_await_error();
3023    else {
3024        GBDATA *gb_data     = GBT_find_sequence(gb_species, test_aliname);
3025        if (!gb_data) error = GB_await_error();
3026        else {
3027            data             = GB_read_char_pntr(gb_data);
3028            if (!data) error = GB_await_error();
3029        }
3030    }
3031
3032    TEST_EXPECT_NULL(error.deliver());
3033   
3034    return data;
3035}
3036
3037static const char *get_used_rels_for(GBDATA *gb_main, const char *species_name) {
3038    GB_transaction  ta(gb_main);
3039    const char     *result     = NULp;
3040    GBDATA         *gb_species = GBT_find_species(gb_main, species_name);
3041    if (!gb_species) result    = GBS_global_string("<No such species '%s'>", species_name);
3042    else {
3043        GBDATA *gb_used_rels      = GB_search(gb_species, "used_rels", GB_FIND);
3044        if (!gb_used_rels) result = "<No such field 'used_rels'>";
3045        else    result            = GB_read_char_pntr(gb_used_rels);
3046    }
3047    return result;
3048}
3049
3050static GB_ERROR forget_used_relatives(GBDATA *gb_main) {
3051    GB_ERROR       error = NULp;
3052    GB_transaction ta(gb_main);
3053    for (GBDATA *gb_species = GBT_first_species(gb_main);
3054         gb_species && !error;
3055         gb_species = GBT_next_species(gb_species))
3056    {
3057        GBDATA *gb_used_rels = GB_search(gb_species, "used_rels", GB_FIND);
3058        if (gb_used_rels) {
3059            error = GB_delete(gb_used_rels);
3060        }
3061    }
3062    return error;
3063}
3064
3065
3066#define ALIGNED_DATA_OF(name) get_aligned_data_of(gb_main, name)
3067#define USED_RELS_FOR(name)   get_used_rels_for(gb_main, name)
3068
3069// ----------------------------------------
3070
3071static GBDATA *selection_fake_gb_main = NULp;
3072static GBDATA *selection_fake_gb_last = NULp;
3073
3074static GBDATA *fake_first_selected(int *count) {
3075    selection_fake_gb_last = NULp;
3076    *count                 = 2; // we fake two species as selected ("c1" and "c2")
3077    return GBT_find_species(selection_fake_gb_main, "c1");
3078}
3079static GBDATA *fake_next_selected() {
3080    if (!selection_fake_gb_last) {
3081        selection_fake_gb_last = GBT_find_species(selection_fake_gb_main, "c2");
3082    }
3083    else {
3084        selection_fake_gb_last = NULp;
3085    }
3086    return selection_fake_gb_last;
3087}
3088
3089static char *fake_get_consensus(const char*, PosRange range) {
3090    const char *data = get_aligned_data_of(selection_fake_gb_main, "s1");
3091    if (range.is_whole()) return strdup(data);
3092    return ARB_strpartdup(data+range.start(), data+range.end());
3093}
3094
3095static void test_install_fakes(GBDATA *gb_main) {
3096    selection_fake_gb_main = gb_main;
3097}
3098
3099
3100// ----------------------------------------
3101
3102static AlignParams test_ali_params = {
3103    FA_NO_REPORT,
3104    false, // showGapsMessages
3105    PosRange::whole()
3106};
3107
3108static AlignParams test_ali_params_partial = {
3109    FA_NO_REPORT,
3110    false, // showGapsMessages
3111    PosRange(25, 60)
3112};
3113
3114// ----------------------------------------
3115
3116static struct arb_unit_test::test_alignment_data TestAlignmentData_TargetAndReferenceHandling[] = {
3117    { 0, "s1", ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........." }, // reference
3118    { 0, "s2", "AUCUCCUAAACCCAACCGUAGUUCGAAUUGAGGACUGUAACUC......................................................" }, // align single sequence (same data as reference)
3119    { 1, "m1", "UAGAGGAUUUGGGUUGGCAUCAAGCUUAACUCCUGACAUUGAG......................................................" }, // align marked sequences.. (complement of reference)
3120    { 1, "m2", "...UCCUAAACCAACCCGUAGUUCGAAUUGAGGACUGUAA........................................................." },
3121    { 1, "m3", "AUC---UAAACCAACCCGUAGUUCGAAUUGAGGACUG---CUC......................................................" },
3122    { 0, "c1", "AUCUCCUAAACCCAACC--------AAUUGAGGACUGUAACUC......................................................" },  // align selected sequences..
3123    { 0, "c2", "AUCUCCU------AACCGUAGUUCCCCGAA------ACUGUAACUC..................................................." },
3124    { 0, "r1", "GAGUUACAGUCCUCAAUUCGGGGAACUACGGUUGGGUUUAGGAGAU..................................................." },  // align by faked pt_server
3125};
3126
3127void TEST_Aligner_TargetAndReferenceHandling() {
3128    // performs some alignments to check logic of target and reference handling
3129   
3130    GB_shell   shell;
3131    ARB_ERROR  error;
3132    GBDATA    *gb_main = TEST_CREATE_DB(error, test_aliname, TestAlignmentData_TargetAndReferenceHandling, false);
3133
3134    TEST_EXPECT_NULL(error.deliver());
3135
3136    SearchRelativeParams search_relative_params(new FakeFamilyFinder(gb_main, test_aliname, false, 8),
3137                                                test_aliname,
3138                                                2);
3139
3140    test_install_fakes(gb_main);
3141    arb_suppress_progress silence;
3142
3143    // bool cont_on_err    = true;
3144    bool cont_on_err = false;
3145
3146    TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 3); // we got 3 marked species
3147    {
3148        Aligner aligner(gb_main,
3149                        FA_CURRENT,
3150                        test_aliname,
3151                        "s2",                       // toalign
3152                        NULp,                       // get_first_selected_species
3153                        NULp,                       // get_next_selected_species
3154                        "s1",                       // reference species
3155                        NULp,                       // get_consensus
3156                        search_relative_params,     // relative search
3157                        FA_TURN_ALWAYS,
3158                        test_ali_params,
3159                        0,
3160                        cont_on_err,
3161                        FA_NO_ACTION);
3162        error = aligner.run();
3163        TEST_EXPECT_NULL(error.deliver());
3164    }
3165    TEST_EXPECT_EQUAL(GBT_count_marked_species(gb_main), 3); // we still got 3 marked species
3166    {
3167        Aligner aligner(gb_main,
3168                        FA_MARKED,
3169                        test_aliname,
3170                        NULp,                       // toalign
3171                        NULp,                       // get_first_selected_species
3172                        NULp,                       // get_next_selected_species
3173                        "s1",                       // reference species
3174                        NULp,                       // get_consensus
3175                        search_relative_params,     // relative search
3176                        FA_TURN_ALWAYS,
3177                        test_ali_params,
3178                        0,
3179                        cont_on_err,
3180                        FA_MARK_FAILED);
3181        error = aligner.run();
3182        TEST_EXPECT_NULL(error.deliver());
3183
3184        TEST_EXPECT(!cont_on_err || GBT_count_marked_species(gb_main) == 0); // FA_MARK_FAILED (none failed -> none marked)
3185    }
3186    {
3187        Aligner aligner(gb_main,
3188                        FA_SELECTED,
3189                        test_aliname,
3190                        NULp,                       // toalign
3191                        fake_first_selected,        // get_first_selected_species
3192                        fake_next_selected,         // get_next_selected_species
3193                        NULp,                       // reference species
3194                        fake_get_consensus,         // get_consensus
3195                        search_relative_params,     // relative search
3196                        FA_TURN_ALWAYS,
3197                        test_ali_params,
3198                        0,
3199                        cont_on_err,
3200                        FA_MARK_ALIGNED);
3201        error = aligner.run();
3202        TEST_EXPECT_NULL(error.deliver());
3203
3204        TEST_EXPECT(!cont_on_err || GBT_count_marked_species(gb_main) == 2); // FA_MARK_ALIGNED (2 selected were aligned)
3205    }
3206    {
3207        Aligner aligner(gb_main,
3208                        FA_CURRENT,
3209                        test_aliname,
3210                        "r1",                       // toalign
3211                        NULp,                       // get_first_selected_species
3212                        NULp,                       // get_next_selected_species
3213                        NULp,                       // reference species
3214                        NULp,                       // get_consensus
3215                        search_relative_params,     // relative search
3216                        FA_TURN_ALWAYS,
3217                        test_ali_params,
3218                        0,
3219                        cont_on_err,
3220                        FA_MARK_ALIGNED);
3221
3222        error = aligner.run();
3223        TEST_EXPECT_NULL(error.deliver());
3224
3225        TEST_EXPECT(!cont_on_err || GBT_count_marked_species(gb_main) == 1);
3226    }
3227   
3228    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C...........");
3229
3230    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".......UAG--AGG-A------U-U-UGGGU-UG-G-C-A-U-CAA-GCU--------UAA-C-UCCUG-AC--A-UUGAG...............");
3231    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), "..............U-C------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA................");
3232    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C...........");
3233
3234    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C...........");
3235    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C-U-AA---------C-C-G-UAG-UUC------------C-CCGAA-AC--U-GUAA-CU-C...........");
3236
3237    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("r1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUCCCC-----GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // here sequence shall be turned!
3238
3239
3240    TEST_EXPECT_EQUAL(USED_RELS_FOR("r1"), "s2:43, s1:1");
3241
3242    // ----------------------------------------------
3243    //      now align all others vs next relative
3244
3245    search_relative_params.maxRelatives = 5;
3246    TEST_EXPECT_NO_ERROR(forget_used_relatives(gb_main));
3247
3248    int species_count = ARRAY_ELEMS(TestAlignmentData_TargetAndReferenceHandling);
3249    for (int sp = 0; sp<species_count; ++sp) {
3250        const char *name = TestAlignmentData_TargetAndReferenceHandling[sp].name;
3251        if (strcmp(name, "r1") != 0) { // skip "r1" (already done above)
3252            Aligner aligner(gb_main,
3253                            FA_CURRENT,
3254                            test_aliname,
3255                            name,                       // toalign
3256                            NULp,                       // get_first_selected_species
3257                            NULp,                       // get_next_selected_species
3258                            NULp,                       // reference species
3259                            NULp,                       // get_consensus
3260                            search_relative_params,     // relative search
3261                            FA_TURN_ALWAYS,
3262                            test_ali_params,
3263                            0,
3264                            cont_on_err,
3265                            FA_MARK_ALIGNED);
3266
3267            error = aligner.run();
3268            TEST_EXPECT_NULL(error.deliver());
3269
3270            TEST_EXPECT(!cont_on_err || GBT_count_marked_species(gb_main) == 1);
3271        }
3272    }
3273
3274    TEST_EXPECT_EQUAL(USED_RELS_FOR("s1"), "s2");
3275    TEST_EXPECT_EQUAL(USED_RELS_FOR("s2"), "s1"); // same as done manually
3276    TEST_EXPECT_EQUAL(USED_RELS_FOR("m1"), "r1:42");
3277    TEST_EXPECT_EQUAL(USED_RELS_FOR("m2"), "m3");
3278    TEST_EXPECT_EQUAL(USED_RELS_FOR("m3"), "m2");
3279    TEST_EXPECT_EQUAL(USED_RELS_FOR("c1"), "r1");
3280    TEST_EXPECT_EQUAL(USED_RELS_FOR("c2"), "r1");
3281
3282    //                                       range aligned below (see test_ali_params_partial)
3283    //                                       "-------------------------XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX------------------------------------"
3284    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // 1st aligning of 's1'
3285    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above (again aligned vs 's1')
3286
3287    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".........U--AGA-G------G---AUUUG-GG-U-U-G-G-CAU-CAAGCU-----UAA-C-UCCUG-AC--A-UUGAG---------------"); // changed; @@@ bug: no dots at end
3288 
3289    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), ".........U--C----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----UA-A..........."); // changed (1st align vs 's1', this align vs 'm3')
3290    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C..........."); // same_as_above (1st align vs 's1', this align vs 'm2')
3291    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3292    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C-U--------A-A-C-C-G-UAG-UUCCCC-----GA--------A-AC--U-GUAA-CU-C..........."); // changed
3293
3294    // --------------------------------------
3295    //      test partial relative search
3296
3297    search_relative_params.getFamilyFinder()->restrict_2_region(test_ali_params_partial.range);
3298    TEST_EXPECT_NO_ERROR(forget_used_relatives(gb_main));
3299
3300    for (int sp = 0; sp<species_count; ++sp) {
3301        const char *name = TestAlignmentData_TargetAndReferenceHandling[sp].name;
3302        Aligner aligner(gb_main,
3303                        FA_CURRENT,
3304                        test_aliname,
3305                        name,                       // toalign
3306                        NULp,                       // get_first_selected_species
3307                        NULp,                       // get_next_selected_species
3308                        NULp,                       // reference species
3309                        NULp,                       // get_consensus
3310                        search_relative_params,     // relative search
3311                        FA_TURN_NEVER,
3312                        test_ali_params_partial,
3313                        0,
3314                        cont_on_err,
3315                        FA_MARK_ALIGNED);
3316
3317        error = aligner.run();
3318        TEST_EXPECT_NULL(error.deliver());
3319
3320        TEST_EXPECT(!cont_on_err || GBT_count_marked_species(gb_main) == 1);
3321    }
3322
3323    TEST_EXPECT_EQUAL(USED_RELS_FOR("s1"), "s2");
3324    TEST_EXPECT_EQUAL(USED_RELS_FOR("s2"), "s1");
3325    TEST_EXPECT_EQUAL(USED_RELS_FOR("m1"), "r1"); // (not really differs)
3326    TEST_EXPECT_EQUAL(USED_RELS_FOR("m2"), "m3");
3327    TEST_EXPECT_EQUAL(USED_RELS_FOR("m3"), "m2");
3328    TEST_EXPECT_EQUAL(USED_RELS_FOR("c1"), "r1");
3329    TEST_EXPECT_EQUAL(USED_RELS_FOR("c2"), "r1");
3330    TEST_EXPECT_EQUAL(USED_RELS_FOR("r1"), "s2:20, c2:3");
3331
3332    //                                       aligned range (see test_ali_params_partial)
3333    //                                       "-------------------------XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX------------------------------------"
3334    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3335    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("s2"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3336
3337    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m1"), ".........U--AGA-G------G-A-UU-UG-GG-U-U-G-G-CAU-CAAGCU-----UAA-C-UCCUG-AC--A-UUGAG---------------"); // changed
3338    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m2"), ".........U--C----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----UA-A..........."); // same_as_above
3339    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("m3"), ".........A--U----------C-U-AAACC-AA-C-C-C-G-UAG-UUC--------GAA-U-UGAGG-AC--U-G----CU-C..........."); // same_as_above
3340
3341    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-------------------AA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3342    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("c2"), ".........A--UCU-C------C---------UA-A-C-C-G-UAG-UUCCCC-----GA--------A-AC--U-GUAA-CU-C..........."); // changed
3343
3344    TEST_EXPECT_EQUAL(ALIGNED_DATA_OF("r1"), ".........A--UCU-C------C-U-AAACC-CA-A-C-C-G-UAG-UUCCCC-----GAA-U-UGAGG-AC--U-GUAA-CU-C..........."); // same_as_above
3345
3346    GB_close(gb_main);
3347}
3348
3349// ----------------------------------------
3350
3351static struct arb_unit_test::test_alignment_data TestAlignmentData_checksumError[] = {
3352    { 0, "MtnK1722", "...G-GGC-C-G............CCC-GG--------CAAUGGGGGCGGCCCGGCGGAC----GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCC....." }, // gets aligned
3353    { 0, "MhnFormi", "...A-CGA-U-C------------CUUCGG--------GGUCG-U-GG-C-GU-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCCU...." }, // 1st relative
3354    { 0, "MhnT1916", "...A-CGA-A-C------------CUU-GU--------GUUCG-U-GG-C-GA-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCCU...." }, // next relative
3355    { 0, "MthVanni", "...U-GGU-U-U------------C-------------GGCCA-U-GG-C-GG-A--C------GG--C-UCAUU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUG-CGGC-UGGAUCACCUCC....." }, // next relative
3356    { 0, "ThcCeler", "...G-GGG-C-G...CC-U---U--------GC--G--CGCAC-C-GG-C-GG-A--C------GG--C-UCAGU-A---AAG-UCGUAACAA-GG-UAG-CCGU-AGGGGAA-CCUA-CGGC-UCGAUCACCUCCU...." }, // next relative
3357};
3358
3359void TEST_SLOW_Aligner_checksumError() {
3360    // @@@ SLOW cause this often gets terminated in nightly builds
3361    //     no idea why. normally it runs a few ms
3362   
3363    // this produces an internal aligner error
3364
3365    GB_shell   shell;
3366    ARB_ERROR  error;
3367    GBDATA    *gb_main = TEST_CREATE_DB(error, test_aliname, TestAlignmentData_checksumError, false);
3368
3369    SearchRelativeParams search_relative_params(new FakeFamilyFinder(gb_main, test_aliname, false, 8),
3370                                                test_aliname,
3371                                                10); // use up to 10 relatives
3372
3373    test_install_fakes(gb_main);
3374    arb_suppress_progress silence;
3375
3376    bool cont_on_err = true;
3377    if (!error) {
3378        Aligner aligner(gb_main,
3379                        FA_CURRENT,
3380                        test_aliname,
3381                        "MtnK1722",                 // toalign
3382                        NULp,                       // get_first_selected_species
3383                        NULp,                       // get_next_selected_species
3384                        NULp,                       // reference species
3385                        NULp,                       // get_consensus
3386                        search_relative_params,     // relative search
3387                        FA_TURN_ALWAYS,
3388                        test_ali_params,
3389                        0,
3390                        cont_on_err,
3391                        FA_MARK_ALIGNED);
3392
3393        error = aligner.run();
3394    }
3395    {
3396        GB_ERROR err = error.deliver();
3397        TEST_EXPECT_NULL__BROKEN(err, "Aligner produced 1 error");
3398    }
3399    TEST_EXPECT_EQUAL__BROKEN(USED_RELS_FOR("MtnK1722"), "???", "<No such field 'used_rels'>"); // subsequent failure
3400
3401    GB_close(gb_main);
3402}
3403
3404static const char *asstr(LooseBases& ub) {
3405    LooseBases tmp;
3406    while (!ub.is_empty()) tmp.memorize(ub.recall());
3407   
3408    const char *result = "";
3409    while (!tmp.is_empty()) {
3410        ExplicitRange r = tmp.recall();
3411        result          = GBS_global_string("%s %i/%i", result, r.start(), r.end());
3412        ub.memorize(r);
3413    }
3414    return result;
3415}
3416
3417void TEST_BASIC_UnalignedBases() {
3418    LooseBases ub;
3419    TEST_EXPECT(ub.is_empty());
3420    TEST_EXPECT_EQUAL(asstr(ub), "");
3421
3422    // test add+remove
3423    ub.memorize(ExplicitRange(5, 7));
3424    TEST_REJECT(ub.is_empty());
3425    TEST_EXPECT_EQUAL(asstr(ub), " 5/7");
3426   
3427    TEST_EXPECT(ub.recall() == ExplicitRange(5, 7));
3428    TEST_EXPECT(ub.is_empty());
3429
3430    ub.memorize(ExplicitRange(2, 4));
3431    TEST_EXPECT_EQUAL(asstr(ub), " 2/4");
3432
3433    ub.memorize(ExplicitRange(4, 9));
3434    TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9");
3435   
3436    ub.memorize(ExplicitRange(8, 10));
3437    ub.memorize(ExplicitRange(11, 14));
3438    ub.memorize(ExplicitRange(12, 17));
3439    TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9 8/10 11/14 12/17");
3440    TEST_EXPECT_EQUAL(asstr(ub), " 2/4 4/9 8/10 11/14 12/17"); // check asstr has no side-effect
3441
3442    {
3443        LooseBases toaddNrecalc;
3444
3445        CompactedSubSequence Old("ACGTACGT", 8, "name1");
3446        CompactedSubSequence New("--A-C--G-T--A-C-G-T", 19, "name2");
3447        // ---------------------- 0123456789012345678
3448
3449        toaddNrecalc.memorize(ExplicitRange(1, 7));
3450        toaddNrecalc.memorize(ExplicitRange(3, 5));
3451        TEST_EXPECT_EQUAL(asstr(toaddNrecalc), " 1/7 3/5");
3452
3453        ub.follow_ali_change_and_append(toaddNrecalc, AliChange(Old, New));
3454
3455        TEST_EXPECT_EQUAL(asstr(ub), " 3/18 8/15 2/4 4/9 8/10 11/14 12/17");
3456        TEST_EXPECT(toaddNrecalc.is_empty());
3457
3458        LooseBases selfRecalc;
3459        selfRecalc.follow_ali_change_and_append(ub, AliChange(New, New));
3460        TEST_EXPECT_EQUAL__BROKEN(asstr(selfRecalc),
3461                                  " 3/18 8/15 0/6 3/11 8/11 10/15 10/17",  // wanted behavior?
3462                                  " 3/18 8/17 0/6 3/11 8/13 10/15 10/18"); // doc wrong behavior @@@ "8/17", "8/13", "10/18" are wrong
3463
3464        ub.follow_ali_change_and_append(selfRecalc, AliChange(New, Old));
3465        TEST_EXPECT_EQUAL__BROKEN(asstr(ub),
3466                                  " 1/7 3/5 0/1 1/3 3/3 4/5 4/6",  // wanted behavior? (from wanted behavior above)
3467                                  " 1/7 3/7 0/2 1/4 3/5 4/6 4/7"); // document wrong behavior
3468        TEST_EXPECT_EQUAL__BROKEN(asstr(ub),
3469                                  " 1/7 3/6 0/1 1/3 3/4 4/5 4/7",  // wanted behavior? (from wrong result above)
3470                                  " 1/7 3/7 0/2 1/4 3/5 4/6 4/7"); // document wrong behavior
3471    }
3472}
3473
3474
3475#endif // UNIT_TESTS
3476
3477
Note: See TracBrowser for help on using the repository browser.