source: tags/testbuild/PARSIMONY/PARS_main.cxx

Last change on this file was 13267, checked in by westram, 11 years ago
  • renamed AP_main::push/pop/clear into remember/revert/accept
  • added convenience wrappers accept_if/revert_if
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 92.6 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : PARS_main.cxx                                     //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "pars_main.hxx"
12#include "pars_klprops.hxx"
13#include "ap_tree_nlen.hxx"
14#include "ap_main.hxx"
15
16#include <aw_awars.hxx>
17#include <aw_preset.hxx>
18#include <aw_msg.hxx>
19#include <aw_root.hxx>
20#include <aw_question.hxx>
21
22#include <awt.hxx>
23#include <awt_sel_boxes.hxx>
24#include <awt_filter.hxx>
25
26#include <ColumnStat.hxx>
27#include <nds.h>
28#include <arb_progress.h>
29#include <arb_misc.h>
30
31#include <gui_aliview.hxx>
32#include <ad_cb.h>
33
34#include <TreeCallbacks.hxx>
35
36#include <list>
37#include <macros.hxx>
38
39#if defined(DEBUG)
40# define TESTMENU
41#endif // DEBUG
42
43using namespace std;
44
45AW_HEADER_MAIN
46
47#define AWAR_COLUMNSTAT_BASE "tmp/pars/colstat"
48#define AWAR_COLUMNSTAT_NAME AWAR_COLUMNSTAT_BASE "/name"
49
50#define AWT_TREE_PARS(ntw) DOWNCAST(AWT_graphic_parsimony*, (ntw)->gfx)
51
52static ArbParsimony *GLOBAL_PARS = NULL;
53
54inline AWT_graphic_parsimony *global_tree() { return GLOBAL_PARS->get_tree(); }
55inline AP_pars_root *global_tree_root() { return global_tree()->get_tree_root(); }
56
57// waaah more globals :(
58AP_main *ap_main; // @@@ move into ArbParsimony? or eliminate ArbParsimony
59
60void ArbParsimony::set_tree(AWT_graphic_parsimony *tree_) {
61    ap_assert(!tree); // only call once
62    tree = tree_;
63    ap_main->set_tree_root(tree);
64}
65
66static void pars_saveNrefresh_changed_tree(AWT_canvas *ntw) {
67    ap_assert((AWT_TREE(ntw) == global_tree()));
68
69    GB_ERROR error = global_tree()->save(ntw->gb_main, 0, 0, 0);
70    if (error) aw_message(error);
71
72    ntw->zoom_reset_and_refresh();
73}
74
75__ATTR__NORETURN static void pars_exit(AW_window *aww) {
76    AW_root *aw_root = aww->get_root();
77    shutdown_macro_recording(aw_root);
78
79    ap_main->accept();
80
81    aw_root->unlink_awars_from_DB(ap_main->get_gb_main());
82#if defined(DEBUG)
83    AWT_browser_forget_db(ap_main->get_gb_main());
84#endif // DEBUG
85    delete ap_main; // closes DB
86    ap_main = NULL;
87
88    exit(EXIT_SUCCESS);
89}
90
91static void AP_user_push_cb(AW_window *aww, AWT_canvas *) {
92    ap_main->user_push();
93    aww->get_root()->awar(AWAR_STACKPOINTER)->write_int(ap_main->get_user_push_counter());
94}
95
96static void AP_user_pop_cb(AW_window *aww, AWT_canvas *ntw) {
97    if (ap_main->get_user_push_counter()<=0) {
98        aw_message("No tree on stack.");
99        return;
100    }
101    ap_main->user_pop();
102    rootNode()->compute_tree();
103    ASSERT_VALID_TREE(rootNode());
104    aww->get_root()->awar(AWAR_STACKPOINTER)->write_int(ap_main->get_user_push_counter());
105
106    pars_saveNrefresh_changed_tree(ntw);
107
108    if (ap_main->get_user_push_counter() <= 0) { // last tree was popped => push again
109        AP_user_push_cb(aww, ntw);
110    }
111}
112
113class InsertData {
114    int  abort_flag;
115    long currentspecies;
116   
117    arb_progress progress;
118
119public:
120
121    bool quick_add_flag;
122    InsertData(bool quick, long spec_count, int add_progress_steps)
123        : abort_flag(false),
124          currentspecies(0),
125          progress(GBS_global_string("Inserting %li species", spec_count), spec_count+add_progress_steps),
126          quick_add_flag(quick)
127    {
128    }
129
130    bool every_sixteenth() { return (currentspecies & 0xf) == 0; }
131
132    bool aborted() const { return abort_flag; }
133    void set_aborted(bool aborted_) { abort_flag = aborted_; }
134
135    void inc() {
136        ++currentspecies;
137        progress.inc();
138        abort_flag = progress.aborted();
139    }
140
141    arb_progress& get_progress() { return progress; }
142};
143
144
145static int sort_sequences_by_length(const char*, long leaf0_ptr, const char*, long leaf1_ptr) {
146    AP_tree *leaf0 = (AP_tree*)leaf0_ptr;
147    AP_tree *leaf1 = (AP_tree*)leaf1_ptr;
148
149    AP_FLOAT len0 = leaf0->get_seq()->weighted_base_count();
150    AP_FLOAT len1 = leaf1->get_seq()->weighted_base_count();
151
152    return len0<len1 ? 1 : (len0>len1 ? -1 : 0); // longest sequence first
153}
154
155static long transform_gbd_to_leaf(const char *key, long val, void *) {
156    if (!val) return val;
157
158#if defined(WARN_TODO)
159#warning use create_linked_leaf() when impl
160#endif
161
162    GBDATA       *gb_node = (GBDATA *)val;
163    AP_pars_root *troot   = ap_main->get_tree_root();
164    AP_tree_nlen *leaf    = DOWNCAST(AP_tree_nlen*, troot->makeNode());
165
166    leaf->forget_origin(); // new leaf is not part of tree yet
167
168    leaf->gb_node = gb_node;
169    leaf->name    = strdup(key);
170    leaf->is_leaf = true;
171
172    leaf->set_seq(troot->get_seqTemplate()->dup());
173    GB_ERROR error = leaf->get_seq()->bind_to_species(gb_node);
174    if (error) {
175        aw_message(error);
176        delete leaf; leaf = 0;
177    }
178    return (long)leaf;
179}
180
181static AP_tree_nlen *insert_species_in_tree(const char *key, AP_tree_nlen *leaf, InsertData *isits) {
182    if (!leaf) return leaf;
183    if (isits->aborted()) return leaf;
184
185    AP_tree_nlen *tree = rootNode();
186
187    if (leaf->get_seq()->weighted_base_count() < MIN_SEQUENCE_LENGTH) {
188        aw_message(GBS_global_string("Species %s has too short sequence (%f, minimum is %i)",
189                                     key,
190                                     leaf->get_seq()->weighted_base_count(),
191                                     MIN_SEQUENCE_LENGTH));
192        delete leaf;
193        return 0;
194    }
195
196    if (!tree) {                                    // no tree yet
197        static AP_tree_nlen *last_inserted = NULL; // @@@ move 'last_inserted' into 'InsertData'
198
199        if (!last_inserted) {                       // store 1st leaf
200            last_inserted = leaf;
201        }
202        else {                                      // 2nd leaf -> create initial tree
203            AP_pars_root *troot = ap_main->get_tree_root();
204
205            leaf->initial_insert(last_inserted, troot);
206            last_inserted = NULL;
207
208            ap_assert(troot->get_root_node() == leaf->get_father());
209            ASSERT_VALID_TREE(troot->get_root_node());
210        }
211    }
212    else {
213        ASSERT_VALID_TREE(tree);
214
215        AP_tree_nlen **branchlist;
216
217        {
218            AP_tree **blist;
219            long      bsum = 0;
220
221            tree->buildBranchList(blist, bsum, true, -1); // get all branches
222            branchlist = (AP_tree_nlen**)blist;
223        }
224
225        AP_tree_nlen *bestposl = tree->get_leftson();
226        AP_tree_nlen *bestposr = tree->get_rightson();
227
228        ap_assert(tree == rootNode());
229        leaf->insert(bestposl);
230        tree = NULL;                                // tree may have changed
231
232        {
233            AP_FLOAT curr_parsimony = rootNode()->costs();
234            AP_FLOAT best_parsimony = curr_parsimony;
235
236            ASSERT_VALID_TREE(rootNode());
237
238            for (long counter = 0; !isits->aborted() && branchlist[counter]; counter += 2) {
239                AP_tree_nlen *bl  = branchlist[counter];
240                AP_tree_nlen *blf = branchlist[counter+1];
241
242                if (blf->father == bl) {
243                    bl  = branchlist[counter+1];
244                    blf = branchlist[counter];
245                }
246
247                ASSERT_VALID_TREE(rootNode());
248                if (bl->father) {
249                    bl->set_root();
250                }
251                ASSERT_VALID_TREE(rootNode());
252
253                leaf->moveNextTo(bl, 0.5);
254                ASSERT_VALID_TREE(rootNode());
255
256                curr_parsimony = rootNode()->costs();
257                ASSERT_VALID_TREE(rootNode());
258
259                if (curr_parsimony < best_parsimony) {
260                    best_parsimony = curr_parsimony;
261                    bestposl = bl;
262                    bestposr = blf;
263                }
264
265            }
266        }
267        delete [] branchlist; branchlist = 0;
268        if (bestposl->father != bestposr) {
269            bestposl = bestposr;
270        }
271
272        ASSERT_VALID_TREE(rootNode());
273
274        leaf->moveNextTo(bestposl, 0.5);
275
276        ASSERT_VALID_TREE(rootNode());
277
278        if (!isits->quick_add_flag) {
279            int deep                           = 5;
280            if (isits->every_sixteenth()) deep = -1;
281           
282            arb_progress progress("optimization");
283            bestposl->get_father()->nn_interchange_rek(deep, AP_BL_NNI_ONLY, true);
284            ASSERT_VALID_TREE(rootNode());
285        }
286        AP_tree_nlen *brother = leaf->get_brother();
287        if (brother->is_leaf && leaf->father->father) {
288            bool brother_is_short = 2 * brother->get_seq()->weighted_base_count() < leaf->get_seq()->weighted_base_count();
289
290            if (brother_is_short) {
291                brother->remove();
292                leaf->remove();
293
294                for (int firstUse = 1; firstUse >= 0; --firstUse) {
295                    AP_tree_nlen *to_insert = firstUse ? leaf : brother;
296                    const char   *format    = firstUse ? "2:%s" : "shortseq:%s";
297
298                    char *label = GBS_global_string_copy(format, to_insert->name);
299                    insert_species_in_tree(label, to_insert, isits);
300                    free(label);
301                }
302            }
303        }
304
305        ASSERT_VALID_TREE(rootNode());
306    }
307
308    return leaf;
309}
310
311static long hash_insert_species_in_tree(const char *key, long leaf, void *cd_isits) {
312    InsertData *isits  = (InsertData*)cd_isits;
313    long        result = (long)insert_species_in_tree(key, (AP_tree_nlen*)leaf, isits);
314    isits->inc();
315    return result;
316}
317
318static long count_hash_elements(const char *, long val, void *cd_count) {
319    if (val) {
320        long *count = (long*)cd_count;
321        (*count)++;
322    }
323    return val;
324}
325
326enum AddWhat {
327    NT_ADD_MARKED,
328    NT_ADD_SELECTED,
329};
330
331static void nt_add(AWT_graphic_parsimony *agt, AddWhat what, bool quick) {
332    GB_ERROR  error = 0;
333
334    AP_tree *oldrootleft  = NULL;
335    AP_tree *oldrootright = NULL;
336    {
337        AP_tree_nlen *root = rootNode();
338        if (root) {
339            root->reset_subtree_layout();
340            oldrootleft  = root->get_leftson();
341            oldrootright = root->get_rightson();
342        }
343    }
344
345    GB_HASH *hash    = 0;
346    GBDATA  *gb_main = agt->gb_main;
347    {
348        GB_transaction ta(gb_main);
349        switch (what) {
350            case NT_ADD_SELECTED: {
351                char *name = GBT_readOrCreate_string(gb_main, AWAR_SPECIES_NAME, "");
352                if (name && strlen(name)) {
353                    GBDATA *gb_species = GBT_find_species(gb_main, name);
354                    if (gb_species) {
355                        hash = GBS_create_hash(1, GB_MIND_CASE);
356                        GBS_write_hash(hash, name, (long)gb_species);
357                    }
358                    else error = GBS_global_string("Selected Species (%s) not found", name);
359                }
360                else error = "Please select a species";
361                free(name);
362                break;
363            }
364            case NT_ADD_MARKED: {
365                hash = GBT_create_marked_species_hash(gb_main);
366                break;
367            }
368        }
369    }
370
371    if (!error) {
372        ap_assert(hash);
373
374        NT_remove_species_in_tree_from_hash(rootNode(), hash);
375
376        long max_species = 0;
377        GBS_hash_do_loop(hash, count_hash_elements, &max_species);
378
379        int        implicitSteps = quick ? 1 : 2; // 1 step for calc_branchlengths, 1 step for NNI
380        InsertData isits(quick, max_species, implicitSteps);
381
382        GB_begin_transaction(gb_main);
383        GBS_hash_do_loop(hash, transform_gbd_to_leaf, NULL);
384        GB_commit_transaction(gb_main);
385
386        {
387            int skipped = max_species - GBS_hash_count_elems(hash);
388            if (skipped) {
389                aw_message(GBS_global_string("Skipped %i species (no data?)", skipped));
390                isits.get_progress().inc_by(skipped);
391            }
392        }
393
394        GBS_hash_do_sorted_loop(hash, hash_insert_species_in_tree, sort_sequences_by_length, &isits);
395
396        if (!quick) {
397            rootEdge()->nni_rek(-1, false, AP_BL_NNI_ONLY, NULL);
398            ++isits.get_progress();
399        }
400
401        if (rootNode()) {
402            rootEdge()->calc_branchlengths();
403            ++isits.get_progress();
404
405            ASSERT_VALID_TREE(rootNode());
406            rootNode()->compute_tree();
407
408            if (oldrootleft) {
409                if (oldrootleft->father == oldrootright) oldrootleft->set_root();
410                else                                     oldrootright->set_root();
411            }
412            else {
413                ARB_edge innermost = rootNode()->get_tree_root()->find_innermost_edge();
414                innermost.set_root();
415            }
416        }
417        else {
418            error = "Tree lost (no leafs left)";
419            isits.get_progress().done();
420        }
421    }
422
423    if (hash) GBS_free_hash(hash);
424    if (error) aw_message(error);
425
426    // @@@ quick-add w/o NNI should sort according to original tree
427    agt->reorder_tree(BIG_BRANCHES_TO_TOP);
428}
429
430// ------------------------------------------
431//      Adding partial sequences to tree
432
433class PartialSequence { 
434    GBDATA               *gb_species;
435    mutable AP_tree_nlen *self;                     // self converted to leaf (ready for insertion)
436    const AP_tree_nlen   *best_full_match;          // full sequence position which matched best
437    long                  overlap;                  // size of overlapping region
438    long                  penalty;                  // weighted mismatches
439    bool                  released;
440    bool                  multi_match;
441    string                multi_list;               // list of equal-rated insertion-points (not containing self)
442
443    AP_tree_nlen *get_self() const {
444        if (!self) {
445            ap_assert(!released); // request not possible, because leaf has already been released!
446
447            self = (AP_tree_nlen*)transform_gbd_to_leaf(GBT_read_name(gb_species), (long)gb_species, NULL);
448            ap_assert(self);
449        }
450        return self;
451    }
452
453public:
454    PartialSequence(GBDATA *gb_species_)
455        : gb_species(gb_species_), self(0), best_full_match(0)
456        , overlap(0),  penalty(LONG_MAX), released(false), multi_match(false)
457    {
458    }
459    PartialSequence(const PartialSequence& other)
460        : gb_species(other.gb_species),
461          self(other.self),
462          best_full_match(other.best_full_match),
463          overlap(other.overlap),
464          penalty(other.penalty),
465          released(other.released),
466          multi_match(other.multi_match),
467          multi_list(other.multi_list)
468    {
469        ap_assert(self == 0); // copying self not implemented
470    }
471    DECLARE_ASSIGNMENT_OPERATOR(PartialSequence);
472    ~PartialSequence() { ap_assert(self == 0); }
473
474    GBDATA *get_species() const { return gb_species; }
475    const AP_tree_nlen *get_best_match() const { return best_full_match; }
476    AP_FLOAT get_branchlength() const { return AP_FLOAT(penalty)/overlap; }
477    void test_match(const AP_tree_nlen *leaf_full);
478    bool is_multi_match() const { return multi_match; }
479
480    const char *get_name() const {
481        const char *name = get_self()->name;
482        ap_assert(name);
483        return name;
484    }
485
486    string get_multilist() const {
487        ap_assert(is_multi_match());
488        return string(best_full_match->name)+multi_list;
489    }
490
491    AP_tree_nlen *release() {
492        AP_tree_nlen *s = self;
493        self            = 0;
494        released        = true;
495        return s;
496    }
497
498    void dump(const char *whichMatch) const {
499        ap_assert(best_full_match);
500        printf("%s match for '%s' is '%s' (overlap=%li penalty=%li)\n",
501               whichMatch, get_name(), best_full_match->name,
502               overlap, penalty);
503    }
504
505};
506
507void PartialSequence::test_match(const AP_tree_nlen *leaf_full) {
508    long curr_overlap;
509    long curr_penalty;
510
511    leaf_full->get_seq()->partial_match(get_self()->get_seq(), &curr_overlap, &curr_penalty);
512
513    bool better = false;
514
515    if (curr_overlap > overlap) {
516        better = true;
517    }
518    else if (curr_overlap == overlap) {
519        if (curr_penalty<penalty) {
520            better = true;
521        }
522        else if (curr_penalty == penalty) {
523            // found two equal-rated insertion points -> store data for warning
524#if defined(DEBUG)
525            if (!multi_match) dump("better");
526            printf("Another equal match is against '%s' (overlap=%li penalty=%li)\n", leaf_full->name, curr_overlap, curr_penalty);
527#endif // DEBUG
528
529            multi_match  = true;
530            multi_list.append(1, '/');
531            multi_list.append(leaf_full->name);
532        }
533    }
534
535    if (better) {
536        overlap         = curr_overlap;
537        penalty         = curr_penalty;
538        best_full_match = leaf_full;
539        multi_match     = false;
540        multi_list      = "";
541
542#if defined(DEBUG)
543        dump("better");
544#endif
545    }
546#if defined(DEBUG)
547    else if (!multi_match) {
548        printf("Worse match against '%s' (overlap=%li penalty=%li)\n", leaf_full->name, curr_overlap, curr_penalty);
549    }
550#endif
551}
552
553static GB_ERROR nt_best_partial_match_rek(list<PartialSequence>& partial, const AP_tree_nlen *tree) {
554    GB_ERROR error = 0;
555
556    if (tree) {
557        if (tree->is_leaf && tree->name) {
558            if (tree->gb_node) {
559                int is_partial = GBT_is_partial(tree->gb_node, 0, true); // marks undef as 'full sequence'
560                if (is_partial == 0) { // do not consider other partial sequences
561                    list<PartialSequence>::iterator i = partial.begin();
562                    list<PartialSequence>::iterator e = partial.end();
563                    for (;  i != e; ++i) {
564                        i->test_match(tree);
565                    }
566                }
567                else if (is_partial == -1) {
568                    error = GB_await_error();
569                }
570            }
571        }
572        else {
573            error             = nt_best_partial_match_rek(partial, tree->get_leftson());
574            if (!error) error = nt_best_partial_match_rek(partial, tree->get_rightson());
575        }
576    }
577    return error;
578}
579
580static void count_partial_and_full(const AP_tree_nlen *at, int *partial, int *full, int *zombies, int default_value, bool define_if_undef) {
581    if (at->is_leaf) {
582        if (at->gb_node) {
583            int is_partial = GBT_is_partial(at->gb_node, default_value, define_if_undef);
584            if (is_partial) ++(*partial);
585            else ++(*full);
586        }
587        else {
588            ++(*zombies);
589        }
590    }
591    else {
592        count_partial_and_full(at->get_leftson(),  partial, full, zombies, default_value, define_if_undef);
593        count_partial_and_full(at->get_rightson(), partial, full, zombies, default_value, define_if_undef);
594    }
595}
596
597static const AP_tree_nlen *find_least_deep_leaf(const AP_tree_nlen *at, int depth, int *min_depth) {
598    if (depth >= *min_depth) {
599        return 0; // already found better or equal
600    }
601
602    if (at->is_leaf) {
603        if (at->gb_node) {
604            *min_depth = depth;
605            return at;
606        }
607        return 0;
608    }
609
610    const AP_tree_nlen *left  = find_least_deep_leaf(at->get_leftson(), depth+1, min_depth);
611    const AP_tree_nlen *right = find_least_deep_leaf(at->get_rightson(), depth+1, min_depth);
612
613    return right ? right : left;
614}
615inline AP_tree_nlen *find_least_deep_leaf(AP_tree_nlen *at, int depth, int *min_depth) {
616    return const_cast<AP_tree_nlen*>(find_least_deep_leaf(const_cast<const AP_tree_nlen*>(at), depth, min_depth));
617}
618
619static long push_partial(const char *, long val, void *cd_partial) {
620    list<PartialSequence> *partial = reinterpret_cast<list<PartialSequence> *>(cd_partial);
621    partial->push_back(PartialSequence((GBDATA*)val));
622    return val;
623}
624
625// -------------------------------
626//      Add Partial sequences
627
628static void nt_add_partial(AWT_graphic_parsimony *agt) {
629    GB_ERROR  error   = NULL;
630    GBDATA   *gb_main = agt->gb_main;
631
632    GB_begin_transaction(gb_main);
633
634    int full_marked_sequences = 0;
635
636    arb_progress part_add_progress("Adding partial sequences");
637
638    {
639        list<PartialSequence> partial;
640        {
641            GB_HASH *partial_hash = GBS_create_hash(GBT_get_species_count(gb_main), GB_MIND_CASE);
642
643            int marked_found             = 0;
644            int partial_marked_sequences = 0;
645            int no_data                  = 0;      // no data in alignment
646
647            for (GBDATA *gb_marked = GBT_first_marked_species(gb_main);
648                 !error && gb_marked;
649                 gb_marked = GBT_next_marked_species(gb_marked))
650            {
651                ++marked_found;
652
653                if (GBT_find_sequence(gb_marked, ap_main->get_aliname())) { // species has sequence in alignment
654                    const char *name = GBT_read_name(gb_marked);
655
656                    switch (GBT_is_partial(gb_marked, 1, true)) { // marks undef as 'partial sequence'
657                        case 0: { // full sequences
658                            GBT_message(gb_main, GBS_global_string("'%s' is a full sequence (cannot add partial)", name));
659                            ++full_marked_sequences;
660                            break;
661                        }
662                        case 1:     // partial sequences
663                            ++partial_marked_sequences;
664                            GBS_write_hash(partial_hash, name, (long)gb_marked);
665                            break;
666                        case -1:    // error
667                            error = GB_await_error();
668                            break;
669                        default:
670                            ap_assert(0);
671                            break;
672                    }
673                }
674                else {
675                    no_data++;
676                }
677            }
678
679            if (!error && !marked_found) error = "There are no marked species";
680
681            if (!error) {
682                NT_remove_species_in_tree_from_hash(rootNode(), partial_hash); // skip all species which are in tree
683                GBS_hash_do_loop(partial_hash, push_partial, &partial); // build partial list from hash
684
685                int partials_already_in_tree = partial_marked_sequences - partial.size();
686
687                if (no_data>0)                  GBT_message(gb_main, GBS_global_string("%i marked species have no data in '%s'",        no_data, ap_main->get_aliname()));
688                if (full_marked_sequences>0)    GBT_message(gb_main, GBS_global_string("%i marked species are declared full sequences", full_marked_sequences));
689                if (partials_already_in_tree>0) GBT_message(gb_main, GBS_global_string("%i marked species are already in tree",         partials_already_in_tree));
690
691                if (partial.empty()) error = "No species left to add";
692            }
693
694            GBS_free_hash(partial_hash);
695        }
696
697        if (!error) error = GBT_add_new_changekey(gb_main, "ARB_partial", GB_INT);
698
699        if (!error) {
700            rootNode()->reset_subtree_layout();
701
702            // find best matching full sequence for each partial sequence
703            error = nt_best_partial_match_rek(partial, rootNode());
704
705            list<PartialSequence>::iterator i = partial.begin();
706            list<PartialSequence>::iterator e = partial.end();
707
708            arb_progress part_insert_progress(partial.size());
709
710#if defined(DEBUG)
711            // show results :
712            for (; i != e; ++i) i->dump("best");
713            i = partial.begin();
714#endif // DEBUG
715
716            for (; i != e && !error; ++i) {
717                const char *name = i->get_name();
718
719                if (i->is_multi_match()) {
720                    GBT_message(gb_main, GBS_global_string("Insertion of '%s' is ambiguous.\n"
721                                                                  "(took first of equal scored insertion points: %s)",
722                                                                  name, i->get_multilist().c_str()));
723                }
724
725                AP_tree_nlen *part_leaf  = i->release();
726                AP_tree_nlen *full_seq   = const_cast<AP_tree_nlen*>(i->get_best_match());
727                AP_tree_nlen *brother    = full_seq->get_brother();
728                int           is_partial = 0;
729                AP_tree_nlen *target     = 0;
730
731                if (brother->is_leaf) {
732                    if (brother->gb_node) {
733                        is_partial = GBT_is_partial(brother->gb_node, 0, true);
734
735                        if (is_partial) { // brother is partial sequence
736                            target = brother; // insert as brother of brother
737                        }
738                        else {
739                            target = full_seq; // insert as brother of full_seq
740                        }
741                    }
742                    else {
743                        error = "There are zombies in your tree - please remove them";
744                    }
745                }
746                else {
747                    int partial_count = 0;
748                    int full_count    = 0;
749                    int zombie_count  = 0;
750
751                    count_partial_and_full(brother, &partial_count, &full_count, &zombie_count, 0, true);
752
753                    if (zombie_count) {
754                        error = "There are zombies in your tree - please remove them";
755                    }
756                    else if (full_count) {
757                        // brother is a subtree containing full sequences
758                        // -> add new brother to full_seq found above
759                        target = full_seq;
760                    }
761                    else {      // brother subtree only contains partial sequences
762                        // find one of the least-deep leafs
763                        int depth  = INT_MAX;
764                        target     = find_least_deep_leaf(brother, 0, &depth);
765                        is_partial = 1;
766                    }
767                }
768
769
770                if (!error) {
771#if defined(DEBUG)
772                    printf("inserting '%s'\n", name);
773#endif // DEBUG
774                    part_leaf->insert(target);
775
776                    // we need to create the sequence of the father node!
777                    AP_tree_nlen *father = part_leaf->get_father();
778                    father->costs();
779
780                    // ensure full-sequence is always on top
781                    if (father->rightson == target) {
782                        father->swap_sons();
783                    }
784
785                    if (!error) { // now correct the branch lengths modified by insert()
786                        // calc the original branchlen (of target leaf branch)
787                        GBT_LEN orglen = father->get_branchlength()+target->get_branchlength();
788
789                        if (is_partial) { // we have a subtree of partial sequences
790                            target->set_branchlength(orglen); // restore original branchlength
791                            father->set_branchlength(0); // all father branches are zero length
792                        }
793                        else { // we have a subtree of one full+one partial sequence
794                            ap_assert(full_seq->get_father() == father);
795
796                            father->set_branchlength(orglen); // father branch represents original length (w/o partial seq)
797                            full_seq->set_branchlength(0);    // full seq has no sub-branch length
798                        }
799                        part_leaf->set_branchlength(i->get_branchlength());
800                        printf("Adding with branchlength=%f\n", i->get_branchlength());
801                    }
802                }
803                else {
804                    delete part_leaf;
805                }
806
807                part_insert_progress.inc_and_check_user_abort(error);
808            }
809        }
810    }
811
812    if (full_marked_sequences) {
813        GBT_message(gb_main, GBS_global_string("%i marked full sequences were not added", full_marked_sequences));
814    }
815
816    if (error) {
817        GBT_message(gb_main, error);
818        GB_abort_transaction(gb_main);
819    }
820    else {
821        GB_commit_transaction(gb_main);
822    }
823}
824
825static void NT_add_partial_and_update(UNFIXED, AWT_canvas *ntw) {
826    nt_add_partial(AWT_TREE_PARS(ntw));
827    pars_saveNrefresh_changed_tree(ntw);
828}
829
830// -------------------------------
831//      add marked / selected
832
833static void nt_add_and_update(AWT_canvas *ntw, AddWhat what, bool quick) {
834    nt_add(AWT_TREE_PARS(ntw), what, quick);
835    pars_saveNrefresh_changed_tree(ntw);
836}
837
838static void NT_add_and_NNI(UNFIXED, AWT_canvas *ntw, AddWhat what) { nt_add_and_update(ntw, what, false); }
839static void NT_add_quick  (UNFIXED, AWT_canvas *ntw, AddWhat what) { nt_add_and_update(ntw, what, true);  }
840
841// ------------------------------------------
842//      remove and add marked / selected
843
844static void nt_reAdd(AWT_graphic_parsimony *agt, AddWhat what, bool quick) {
845    if (agt->get_root_node()) {
846        ap_assert(what == NT_ADD_MARKED); // code below will misbehave for NT_ADD_SELECTED
847        agt->get_tree_root()->remove_leafs(AWT_RemoveType(AWT_REMOVE_BUT_DONT_FREE|AWT_REMOVE_MARKED));
848        nt_add(agt, what, quick);
849    }
850}
851
852static void nt_reAdd_and_update(AWT_canvas *ntw, AddWhat what, bool quick) {
853    nt_reAdd(AWT_TREE_PARS(ntw), what, quick);
854    pars_saveNrefresh_changed_tree(ntw);
855}
856
857static void NT_reAdd_and_NNI(UNFIXED, AWT_canvas *ntw, AddWhat what) { nt_reAdd_and_update(ntw, what, false); }
858static void NT_reAdd_quick  (UNFIXED, AWT_canvas *ntw, AddWhat what) { nt_reAdd_and_update(ntw, what, true);  }
859
860// --------------------------------------------------------------------------------
861
862static void calc_branchlengths(AWT_graphic_parsimony *agt) {
863    arb_progress progress("Calculating Branch Lengths");
864    rootEdge()->calc_branchlengths();
865    agt->reorder_tree(BIG_BRANCHES_TO_TOP);
866}
867
868static void NT_calc_branch_lengths(AW_window *, AWT_canvas *ntw) {
869    calc_branchlengths(AWT_TREE_PARS(ntw));
870    pars_saveNrefresh_changed_tree(ntw);
871}
872
873static void NT_bootstrap(AW_window *, AWT_canvas *ntw, bool limit_only) {
874    arb_progress progress("Calculating Bootstrap Limit");
875    AP_BL_MODE mode       = AP_BL_MODE((limit_only ? AP_BL_BOOTSTRAP_LIMIT : AP_BL_BOOTSTRAP_ESTIMATE)|AP_BL_BL_ONLY);
876    rootEdge()->nni_rek(-1, false, mode, NULL);
877    AWT_TREE(ntw)->reorder_tree(BIG_BRANCHES_TO_TOP);
878    AWT_TREE(ntw)->displayed_root = AWT_TREE(ntw)->get_root_node();
879    pars_saveNrefresh_changed_tree(ntw);
880}
881
882static void optimizeTree(AWT_graphic_parsimony *agt) {
883    arb_progress progress("Optimizing Tree");
884    agt->get_parsimony().optimize_tree(rootNode(), progress);
885    ASSERT_VALID_TREE(rootNode());
886    rootEdge()->calc_branchlengths();
887    agt->reorder_tree(BIG_BRANCHES_TO_TOP);
888    rootNode()->compute_tree();
889}
890static void NT_optimize(AW_window *, AWT_canvas *ntw) {
891    optimizeTree(AWT_TREE_PARS(ntw));
892    pars_saveNrefresh_changed_tree(ntw);
893}
894
895static void recursiveNNI(AWT_graphic_parsimony *agt) {
896    arb_progress progress("Recursive NNI");
897    AP_FLOAT orgPars = rootNode()->costs();
898    AP_FLOAT prevPars = orgPars;
899    progress.subtitle(GBS_global_string("Old parsimony: %.1f", orgPars));
900    while (!progress.aborted()) {
901        AP_FLOAT currPars = rootEdge()->nni_rek(-1, true, AP_BL_NNI_ONLY, NULL);
902        if (currPars == prevPars) break; // no improvement -> abort
903        progress.subtitle(GBS_global_string("New parsimony: %.1f (gain: %.1f)", currPars, orgPars-currPars));
904        prevPars = currPars;
905    }
906    rootEdge()->calc_branchlengths();
907    agt->reorder_tree(BIG_BRANCHES_TO_TOP);
908    rootNode()->compute_tree();
909}
910
911static void NT_recursiveNNI(AW_window *, AWT_canvas *ntw) {
912    recursiveNNI(AWT_TREE_PARS(ntw));
913    pars_saveNrefresh_changed_tree(ntw);
914}
915
916static AW_window *PARS_create_tree_settings_window(AW_root *aw_root) {
917    static AW_window_simple *aws = 0;
918    if (!aws) {
919        aws = new AW_window_simple;
920        aws->init(aw_root, "SAVE_DB", "SAVE ARB DB");
921        aws->load_xfig("awt/tree_settings.fig");
922
923        aws->at("close"); aws->callback((AW_CB0)AW_POPDOWN);
924        aws->create_button("CLOSE", "CLOSE", "C");
925
926        aws->at("help"); aws->callback(makeHelpCallback("nt_tree_settings.hlp"));
927        aws->create_button("HELP", "HELP", "H");
928
929        aws->at("button");
930        aws->auto_space(10, 10);
931        aws->label_length(30);
932
933        aws->label("Base Line Width");
934        aws->create_input_field(AWAR_DTREE_BASELINEWIDTH, 4);
935        aws->at_newline();
936
937        aws->label("Relative vert. Dist");
938        aws->create_input_field(AWAR_DTREE_VERICAL_DIST, 4);
939        aws->at_newline();
940
941        TREE_insert_jump_option_menu(aws, "On species change", AWAR_DTREE_AUTO_JUMP);
942        TREE_insert_jump_option_menu(aws, "On tree change",    AWAR_DTREE_AUTO_JUMP_TREE);
943    }
944    return aws;
945}
946
947// -----------------------
948//      test functions
949
950#if defined(TESTMENU)
951static void refreshTree(AWT_canvas *ntw) {
952    GB_transaction ta(ntw->gb_main);
953
954    AWT_TREE(ntw)->check_update(ntw->gb_main);
955    GB_ERROR error = AWT_TREE(ntw)->save(ntw->gb_main, 0, 0, 0);
956    if (error) aw_message(error);
957    ntw->zoom_reset_and_refresh();
958}
959#endif // TESTMENU
960
961static void setBranchlens(AP_tree_nlen *node, double newLen)
962{
963    node->setBranchlen(newLen, newLen);
964
965    if (!node->is_leaf)
966    {
967        setBranchlens(node->get_leftson(), newLen);
968        setBranchlens(node->get_rightson(), newLen);
969    }
970}
971#if defined(TESTMENU)
972static void TESTMENU_setBranchlen(AW_window *, AWT_canvas *ntw)
973{
974    AP_tree_nlen *root = rootNode();
975
976    setBranchlens(root, 1.0);
977    refreshTree(ntw);
978}
979
980static void TESTMENU_treeStats(AW_window *) {
981    ARB_tree_info tinfo;
982    AP_tree_nlen *root = rootNode();
983
984    if (root) {
985        {
986            GB_transaction ta(root->get_tree_root()->get_gb_main());
987            root->calcTreeInfo(tinfo);
988        }
989
990        puts("Tree stats:");
991
992        printf("nodes      =%6zu\n", tinfo.nodes());
993        printf(" inner     =%6zu\n", tinfo.innerNodes);
994        printf("  groups   =%6zu\n", tinfo.groups);
995        printf(" leafs     =%6zu\n", tinfo.leafs);
996        printf("  unlinked =%6zu (zombies?)\n", tinfo.unlinked);
997        printf("  linked   =%6zu\n", tinfo.linked());
998        printf("   marked  =%6zu\n", tinfo.marked);
999    }
1000    else {
1001        puts("No tree");
1002    }
1003}
1004
1005static void TESTMENU_mixTree(AW_window *, AWT_canvas *ntw)
1006{
1007    rootEdge()->mixTree(100);
1008    refreshTree(ntw);
1009}
1010
1011static void TESTMENU_sortTreeByName(AW_window *, AWT_canvas *ntw)
1012{
1013    AP_tree_nlen *root = rootNode();
1014
1015    root->sortByName();
1016    refreshTree(ntw);
1017}
1018
1019static void TESTMENU_buildAndDumpChain(AW_window *)
1020{
1021    AP_tree_nlen *root = rootNode();
1022
1023    root->get_leftson()->edgeTo(root->get_rightson())->testChain(2);
1024    root->get_leftson()->edgeTo(root->get_rightson())->testChain(3);
1025    root->get_leftson()->edgeTo(root->get_rightson())->testChain(4);
1026    root->get_leftson()->edgeTo(root->get_rightson())->testChain(5);
1027    root->get_leftson()->edgeTo(root->get_rightson())->testChain(6);
1028    root->get_leftson()->edgeTo(root->get_rightson())->testChain(7);
1029    root->get_leftson()->edgeTo(root->get_rightson())->testChain(8);
1030    root->get_leftson()->edgeTo(root->get_rightson())->testChain(9);
1031    root->get_leftson()->edgeTo(root->get_rightson())->testChain(10);
1032    root->get_leftson()->edgeTo(root->get_rightson())->testChain(11);
1033    root->get_leftson()->edgeTo(root->get_rightson())->testChain(-1);
1034}
1035
1036static void init_TEST_menu(AW_window_menu_modes *awm, AWT_canvas *ntw)
1037{
1038    awm->create_menu("Test[debug]", "g", AWM_ALL);
1039
1040    awm->insert_menu_topic("mixtree",         "Mix tree",           "M", "", AWM_ALL, makeWindowCallback(TESTMENU_mixTree, ntw));
1041    awm->insert_menu_topic("treestat",        "Tree statistics",    "s", "", AWM_ALL, TESTMENU_treeStats);
1042    awm->insert_menu_topic("setlens",         "Set branchlens",     "b", "", AWM_ALL, makeWindowCallback(TESTMENU_setBranchlen, ntw));
1043    awm->insert_menu_topic("sorttreebyname",  "Sort tree by name",  "o", "", AWM_ALL, makeWindowCallback(TESTMENU_sortTreeByName, ntw));
1044    awm->insert_menu_topic("buildndumpchain", "Build & dump chain", "c", "", AWM_ALL, TESTMENU_buildAndDumpChain);
1045}
1046#endif // TESTMENU
1047
1048static GB_ERROR pars_check_size(AW_root *awr, GB_ERROR& warning) {
1049    GB_ERROR error = NULL;
1050    warning        = NULL;
1051
1052    char *tree_name = awr->awar(AWAR_TREE)->read_string();
1053    char *filter    = awr->awar(AWAR_FILTER_FILTER)->read_string();
1054    long  ali_len   = 0;
1055
1056    if (strlen(filter)) {
1057        int i;
1058        for (i=0; filter[i]; i++) {
1059            if (filter[i] != '0') ali_len++;
1060        }
1061    }
1062    else {
1063        char *ali_name = awr->awar(AWAR_ALIGNMENT)->read_string();
1064        ali_len        = GBT_get_alignment_len(ap_main->get_gb_main(), ali_name);
1065        if (ali_len<=0) {
1066            error = "Please select a valid alignment";
1067            GB_clear_error();
1068        }
1069        free(ali_name);
1070    }
1071
1072    if (!error) {
1073        long tree_size = GBT_size_of_tree(ap_main->get_gb_main(), tree_name);
1074        if (tree_size == -1) {
1075            error = "Please select an existing tree";
1076        }
1077        else {
1078            unsigned long expected_memuse = (ali_len * tree_size * 4 / 1024);
1079            if (expected_memuse > GB_get_usable_memory()) {
1080                warning = GBS_global_string("Estimated memory usage (%s) exceeds physical memory (will swap)\n"
1081                                            "(did you specify a filter?)",
1082                                            GBS_readable_size(expected_memuse, "b"));
1083            }
1084        }
1085    }
1086
1087    free(filter);
1088    free(tree_name);
1089
1090    ap_assert(!GB_have_error());
1091    return error;
1092}
1093
1094static void pars_reset_optimal_parsimony(AW_window *aww, AWT_canvas *ntw) {
1095    AW_root *awr = aww->get_root();
1096    awr->awar(AWAR_BEST_PARSIMONY)->write_int(awr->awar(AWAR_PARSIMONY)->read_int());
1097    ntw->refresh();
1098}
1099
1100
1101static void pars_start_cb(AW_window *aw_parent, WeightedFilter *wfilt, const PARS_commands *cmds) {
1102    AW_root *awr     = aw_parent->get_root();
1103    GBDATA  *gb_main = ap_main->get_gb_main();
1104    GB_begin_transaction(gb_main);
1105    {
1106        GB_ERROR warning;
1107        GB_ERROR error = pars_check_size(awr, warning);
1108
1109        if (warning && !error) {
1110            char *question = GBS_global_string_copy("%s\nDo you want to continue?", warning);
1111            bool  cont     = aw_ask_sure("swap_warning", question);
1112            free(question);
1113
1114            if (!cont) error = "User abort";
1115
1116        }
1117
1118        if (error) {
1119            aw_message(error);
1120            GB_commit_transaction(gb_main);
1121            return;
1122        }
1123    }
1124
1125
1126    AW_window_menu_modes *awm = new AW_window_menu_modes;
1127    awm->init(awr, "ARB_PARSIMONY", "ARB_PARSIMONY", 400, 200);
1128
1129    GLOBAL_PARS->generate_tree(wfilt);
1130
1131    AWT_canvas *ntw;
1132    {
1133        AP_tree_display_type  old_sort_type = global_tree()->tree_sort;
1134        global_tree()->set_tree_type(AP_LIST_SIMPLE, NULL); // avoid NDS warnings during startup
1135        ntw = new AWT_canvas(gb_main, awm, awm->get_window_id(), global_tree(), AWAR_TREE);
1136        global_tree()->set_tree_type(old_sort_type, ntw);
1137    }
1138
1139    {
1140        GB_ERROR error = 0;
1141        arb_progress progress("loading tree");
1142        NT_reload_tree_event(awr, ntw, 0);             // load tree (but do not expose - first zombies need to be removed)
1143        if (!global_tree()->get_root_node()) {
1144            error = "Failed to load the selected tree";
1145        }
1146        else {
1147            AP_tree_edge::initialize(rootNode());   // builds edges
1148            long removed = global_tree_root()->remove_leafs(AWT_REMOVE_ZOMBIES);
1149
1150            PARS_tree_init(global_tree());
1151            removed += global_tree_root()->remove_leafs(AWT_RemoveType(AWT_REMOVE_ZOMBIES | AWT_REMOVE_NO_SEQUENCE));
1152
1153            if (!global_tree()->get_root_node()) {
1154                const char *aliname = global_tree_root()->get_aliview()->get_aliname();
1155                error               = GBS_global_string("Less than 2 species contain data in '%s'\n"
1156                                                        "Tree vanished", aliname);
1157            }
1158            else if (removed) {
1159                aw_message(GBS_global_string("Removed %li leafs (zombies or species w/o data in alignment)", removed));
1160            }
1161
1162            error = GB_end_transaction(ntw->gb_main, error);
1163            if (!error) {
1164                progress.subtitle("Calculating inner nodes");
1165                GLOBAL_PARS->get_root_node()->costs();
1166            }
1167        }
1168        if (error) aw_popup_exit(error);
1169    }
1170
1171    awr->awar(AWAR_COLOR_GROUPS_USE)->add_callback(makeRootCallback(TREE_recompute_cb, ntw));
1172
1173    if (cmds->add_marked)           NT_add_quick(NULL, ntw, NT_ADD_MARKED);
1174    if (cmds->add_selected)         NT_add_quick(NULL, ntw, NT_ADD_SELECTED);
1175    if (cmds->calc_branch_lengths)  NT_calc_branch_lengths(awm, ntw);
1176    if (cmds->calc_bootstrap)       NT_bootstrap(awm, ntw, 0);
1177    if (cmds->quit)                 pars_exit(awm);
1178
1179    GB_transaction ta(ntw->gb_main);
1180
1181    GBDATA *gb_arb_presets = GB_search(ntw->gb_main, "arb_presets", GB_CREATE_CONTAINER);
1182    GB_add_callback(gb_arb_presets, GB_CB_CHANGED, makeDatabaseCallback(AWT_expose_cb, ntw));
1183
1184#if defined(DEBUG)
1185    AWT_create_debug_menu(awm);
1186#endif // DEBUG
1187
1188    awm->create_menu("File", "F", AWM_ALL);
1189    {
1190        insert_macro_menu_entry(awm, false);
1191        awm->insert_menu_topic("print_tree", "Print Tree ...", "P", "tree2prt.hlp", AWM_ALL, makeWindowCallback(AWT_popup_print_window, ntw));
1192        awm->insert_menu_topic("quit",       "Quit",           "Q", "quit.hlp",     AWM_ALL, pars_exit);
1193    }
1194
1195    awm->create_menu("Species", "S", AWM_ALL);
1196    {
1197        NT_insert_mark_submenus(awm, ntw, 0);
1198
1199    }
1200    awm->create_menu("Tree", "T", AWM_ALL);
1201    {
1202
1203        awm->insert_menu_topic("nds",       "NDS (Node Display Setup) ...",      "N", "props_nds.hlp",   AWM_ALL, makeCreateWindowCallback(AWT_create_nds_window, ntw->gb_main));
1204
1205        awm->sep______________();
1206        awm->insert_menu_topic("tree_2_xfig", "Export tree to XFIG ...", "E", "tree2file.hlp", AWM_ALL, makeWindowCallback(AWT_popup_tree_export_window, ntw));
1207        awm->insert_menu_topic("tree_print",  "Print tree ...",          "P", "tree2prt.hlp",  AWM_ALL, makeWindowCallback(AWT_popup_print_window,       ntw));
1208        awm->sep______________();
1209        NT_insert_collapse_submenu(awm, ntw);
1210        awm->sep______________();
1211        awm->insert_sub_menu("Remove Species from Tree",     "R");
1212        {
1213            awm->insert_menu_topic("tree_remove_deleted", "Remove Zombies", "Z", "trm_del.hlp",    AWM_ALL, makeWindowCallback(NT_remove_leafs, ntw, AWT_RemoveType(AWT_REMOVE_BUT_DONT_FREE|AWT_REMOVE_ZOMBIES)));
1214            awm->insert_menu_topic("tree_remove_marked",  "Remove Marked",  "M", "trm_mrkd.hlp",   AWM_ALL, makeWindowCallback(NT_remove_leafs, ntw, AWT_RemoveType(AWT_REMOVE_BUT_DONT_FREE|AWT_REMOVE_MARKED)));
1215            awm->insert_menu_topic("tree_keep_marked",    "Keep Marked",    "K", "tkeep_mrkd.hlp", AWM_ALL, makeWindowCallback(NT_remove_leafs, ntw, AWT_RemoveType(AWT_REMOVE_BUT_DONT_FREE|AWT_KEEP_MARKED)));
1216        }
1217        awm->close_sub_menu();
1218        awm->insert_sub_menu("Add Species to Tree",      "A");
1219        {
1220            awm->insert_menu_topic("add_marked",         "Add Marked Species",                              "M", "pa_quick.hlp",   AWM_ALL, makeWindowCallback(NT_add_quick,     ntw, NT_ADD_MARKED));
1221            awm->insert_menu_topic("add_marked_nni",     "Add Marked Species + Local Optimization (NNI)",   "N", "pa_add.hlp",     AWM_ALL, makeWindowCallback(NT_add_and_NNI,   ntw, NT_ADD_MARKED));
1222            awm->insert_menu_topic("rm_add_marked",      "Remove & Add Marked Species",                     "R", "pa_quick.hlp",   AWM_ALL, makeWindowCallback(NT_reAdd_quick,   ntw, NT_ADD_MARKED));
1223            awm->insert_menu_topic("rm_add_marked_nni|", "Remove & Add Marked + Local Optimization (NNI)",  "L", "pa_add.hlp",     AWM_ALL, makeWindowCallback(NT_reAdd_and_NNI, ntw, NT_ADD_MARKED));
1224            awm->sep______________();
1225            awm->insert_menu_topic("add_marked_partial", "Add Marked Partial Species",                      "P", "pa_partial.hlp", AWM_ALL, makeWindowCallback(NT_add_partial_and_update, ntw));
1226            awm->sep______________();
1227            awm->insert_menu_topic("add_selected",       "Add Selected Species",                            "S", "pa_quick.hlp",   AWM_ALL, makeWindowCallback(NT_add_quick,     ntw, NT_ADD_SELECTED));
1228            awm->insert_menu_topic("add_selected_nni",   "Add Selected Species + Local Optimization (NNI)", "O", "pa_add.hlp",     AWM_ALL, makeWindowCallback(NT_add_and_NNI,   ntw, NT_ADD_SELECTED));
1229        }
1230        awm->close_sub_menu();
1231        awm->sep______________();
1232        awm->insert_sub_menu("Tree Optimization",        "O");
1233        {
1234            awm->insert_menu_topic("nni",             "Local Optimization (NNI) of Marked Visible Nodes", "L", "",                 AWM_ALL, makeWindowCallback(NT_recursiveNNI, ntw));
1235            awm->insert_menu_topic("kl_optimization", "Global Optimization of Marked Visible Nodes",      "G", "pa_optimizer.hlp", AWM_ALL, makeWindowCallback(NT_optimize,     ntw));
1236        }
1237        awm->close_sub_menu();
1238        awm->insert_menu_topic("reset", "Reset optimal parsimony", "s", "", AWM_ALL, makeWindowCallback(pars_reset_optimal_parsimony, ntw));
1239        awm->sep______________();
1240        awm->insert_menu_topic("beautify_tree",       "Beautify Tree",            "B", "resorttree.hlp",       AWM_ALL, makeWindowCallback(NT_resort_tree_cb, ntw, BIG_BRANCHES_TO_TOP));
1241        awm->insert_menu_topic("calc_branch_lengths", "Calculate Branch Lengths", "L", "pa_branchlengths.hlp", AWM_ALL, makeWindowCallback(NT_calc_branch_lengths, ntw));
1242        awm->sep______________();
1243        awm->insert_menu_topic("calc_upper_bootstrap_indep", "Calculate Upper Bootstrap Limit (dependent NNI)",   "d", "pa_bootstrap.hlp", AWM_ALL, makeWindowCallback(NT_bootstrap,        ntw, false));
1244        awm->insert_menu_topic("calc_upper_bootstrap_dep",   "Calculate Upper Bootstrap Limit (independent NNI)", "i", "pa_bootstrap.hlp", AWM_ALL, makeWindowCallback(NT_bootstrap,        ntw, true));
1245        awm->insert_menu_topic("tree_remove_remark",         "Remove Bootstrap Values",                           "V", "trm_boot.hlp",     AWM_ALL, makeWindowCallback(NT_remove_bootstrap, ntw));
1246    }
1247
1248#if defined(TESTMENU)
1249    init_TEST_menu(awm, ntw);
1250#endif // TESTMENU
1251
1252    awm->create_menu("Reset", "R", AWM_ALL);
1253    {
1254        awm->insert_menu_topic("reset_logical_zoom",  "Logical Zoom",  "L", "rst_log_zoom.hlp",  AWM_ALL, makeWindowCallback(NT_reset_lzoom_cb, ntw));
1255        awm->insert_menu_topic("reset_physical_zoom", "Physical Zoom", "P", "rst_phys_zoom.hlp", AWM_ALL, makeWindowCallback(NT_reset_pzoom_cb, ntw));
1256    }
1257
1258    awm->create_menu("Properties", "P", AWM_ALL);
1259    {
1260        awm->insert_menu_topic("props_menu",  "Menu: Colors and Fonts ...", "M", "props_frame.hlp",      AWM_ALL, AW_preset_window);
1261        awm->insert_menu_topic("props_tree",  "Tree: Colors and Fonts ...", "C", "pars_props_data.hlp",  AWM_ALL, makeCreateWindowCallback(AW_create_gc_window, ntw->gc_manager));
1262        awm->insert_menu_topic("props_tree2", "Tree: Settings ...",         "T", "nt_tree_settings.hlp", AWM_ALL, PARS_create_tree_settings_window);
1263        awm->insert_menu_topic("props_kl",    "KERN. LIN ...",              "K", "kernlin.hlp",          AWM_ALL, makeCreateWindowCallback(create_kernighan_window));
1264        awm->sep______________();
1265        AW_insert_common_property_menu_entries(awm);
1266        awm->sep______________();
1267        awm->insert_menu_topic("save_props", "Save Defaults (pars.arb)", "D", "savedef.hlp", AWM_ALL, AW_save_properties);
1268    }
1269    awm->button_length(5);
1270
1271    awm->insert_help_topic("ARB_PARSIMONY help", "P", "arb_pars.hlp", AWM_ALL, makeHelpCallback("arb_pars.hlp"));
1272
1273    // ----------------------
1274    //      mode buttons
1275    //
1276    // keep them synchronized as far as possible with those in ARB_PARSIMONY
1277    // see ../NTREE/NT_extern.cxx@keepModesSynchronized
1278
1279    awm->create_mode("mode_select.xpm", "mode_select.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_SELECT));
1280    awm->create_mode("mode_mark.xpm",   "mode_mark.hlp",   AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_MARK));
1281    awm->create_mode("mode_group.xpm",  "mode_group.hlp",  AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_GROUP));
1282    awm->create_mode("mode_zoom.xpm",   "mode_pzoom.hlp",  AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_ZOOM));
1283    awm->create_mode("mode_lzoom.xpm",  "mode_lzoom.hlp",  AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_LZOOM));
1284
1285    // reserve mode-locations (to put the modes below at the same position as in ARB_NT)
1286    awm->create_mode("mode_empty.xpm", "mode.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_EMPTY));
1287    awm->create_mode("mode_empty.xpm", "mode.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_EMPTY));
1288
1289    // topology-modification-modes
1290    awm->create_mode("mode_setroot.xpm", "mode_setroot.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_SETROOT));
1291    awm->create_mode("mode_swap.xpm",    "mode_swap.hlp",    AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_SWAP));
1292    awm->create_mode("mode_move.xpm",    "mode_move.hlp",    AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_MOVE));
1293
1294    awm->create_mode("mode_nni.xpm",      "mode_nni.hlp",      AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_NNI));
1295    awm->create_mode("mode_kernlin.xpm",  "mode_kernlin.hlp",  AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_KERNINGHAN));
1296    awm->create_mode("mode_optimize.xpm", "mode_optimize.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_OPTIMIZE));
1297
1298    awm->at(5, 2);
1299    awm->auto_space(0, -2);
1300    awm->shadow_width(1);
1301
1302
1303    int db_treex, db_treey;
1304    awm->get_at_position(&db_treex, &db_treey);
1305    awm->callback(makeHelpCallback("nt_tree_select.hlp"));
1306    awm->button_length(16);
1307    awm->help_text("nt_tree_select.hlp");
1308    awm->create_button(0, AWAR_TREE);
1309
1310
1311    int db_stackx, db_stacky;
1312    awm->label_length(8);
1313    awm->label("Stored");
1314    awm->get_at_position(&db_stackx, &db_stacky);
1315    awm->button_length(6);
1316    awm->callback(makeHelpCallback("ap_stack.hlp"));
1317    awm->help_text("ap_stack.hlp");
1318    awm->create_button(0, AWAR_STACKPOINTER);
1319
1320    int db_parsx, db_parsy;
1321    awm->label_length(14);
1322    awm->label("Current Pars:");
1323    awm->get_at_position(&db_parsx, &db_parsy);
1324
1325    awm->button_length(10);
1326    awm->create_button(0, AWAR_PARSIMONY, 0, "+");
1327
1328    awm->button_length(0);
1329
1330    awm->callback(makeWindowCallback(NT_jump_cb, ntw, AP_JUMP_BY_BUTTON));
1331    awm->help_text("tr_jump.hlp");
1332    awm->create_button("JUMP", "Jump", 0);
1333
1334    awm->callback(makeHelpCallback("arb_pars.hlp"));
1335    awm->help_text("help.hlp");
1336    awm->create_button("HELP", "HELP", "H");
1337
1338    awm->at_newline();
1339
1340    awm->button_length(8);
1341    awm->at_x(db_stackx);
1342    awm->callback(makeWindowCallback(AP_user_pop_cb, ntw));
1343    awm->help_text("ap_stack.hlp");
1344    awm->create_button("POP", "RESTORE", 0);
1345
1346    awm->button_length(6);
1347    awm->callback(makeWindowCallback(AP_user_push_cb, ntw));
1348    awm->help_text("ap_stack.hlp");
1349    awm->create_button("PUSH", "STORE", 0);
1350
1351    awm->at_x(db_parsx);
1352    awm->label_length(14);
1353    awm->label("Optimal Pars:");
1354
1355    awm->button_length(10);
1356    awm->create_button(0, AWAR_BEST_PARSIMONY, 0, "+");
1357
1358    awm->button_length(0);
1359    awm->auto_space(0, -2);
1360
1361    awm->at_x(db_treex);
1362    awm->callback(makeWindowCallback(NT_set_tree_style, ntw, AP_TREE_RADIAL));
1363    awm->help_text("tr_type_radial.hlp");
1364    awm->create_button("RADIAL_TREE", "#radial.xpm", 0);
1365
1366    awm->callback(makeWindowCallback(NT_set_tree_style, ntw, AP_TREE_NORMAL));
1367    awm->help_text("tr_type_list.hlp");
1368    awm->create_button("LIST_TREE", "#dendro.xpm", 0);
1369
1370    awm->at_newline();
1371    awm->at(db_treex, awm->get_at_yposition());
1372
1373    {
1374        SmartPtr<AW_at_storage> maxSize(AW_at_storage::make(awm, AW_AT_MAXSIZE));
1375
1376        awm->button_length(AWAR_FOOTER_MAX_LEN);
1377        awm->create_button(0, AWAR_FOOTER);
1378        awm->at_newline();
1379        awm->restore_at_from(*maxSize);
1380    }
1381
1382    awm->get_at_position(&db_treex, &db_treey);
1383    awm->set_info_area_height(db_treey);
1384
1385    awm->set_bottom_area_height(0);
1386
1387    aw_parent->hide(); // hide parent
1388    awm->show();
1389
1390    awr->awar(AWAR_SPECIES_NAME)->add_callback(makeRootCallback(TREE_auto_jump_cb, ntw, false));
1391
1392    AP_user_push_cb(aw_parent, ntw); // push initial tree
1393}
1394
1395static AW_window *create_pars_init_window(AW_root *awr, const PARS_commands *cmds) {
1396    AW_window_simple *aws = new AW_window_simple;
1397    aws->init(awr, "PARS_PROPS", "SET PARSIMONY OPTIONS");
1398    aws->load_xfig("pars/init.fig");
1399
1400    aws->button_length(10);
1401    aws->label_length(10);
1402
1403    aws->callback(pars_exit);
1404    aws->at("close");
1405    aws->create_button("ABORT", "ABORT", "A");
1406
1407    aws->callback(makeHelpCallback("arb_pars_init.hlp"));
1408    aws->at("help");
1409    aws->create_button("HELP", "HELP", "H");
1410
1411    GBDATA *gb_main                 = ap_main->get_gb_main();
1412    WeightedFilter *weighted_filter = // do NOT free (bound to callbacks)
1413        new WeightedFilter(gb_main, aws->get_root(), AWAR_FILTER_NAME, AWAR_COLUMNSTAT_NAME, aws->get_root()->awar_string(AWAR_ALIGNMENT));
1414
1415    aws->at("filter");
1416    aws->callback(makeCreateWindowCallback(awt_create_select_filter_win, weighted_filter->get_adfiltercbstruct()));
1417    aws->create_button("SELECT_FILTER", AWAR_FILTER_NAME);
1418
1419    aws->at("weights");
1420    aws->callback(makeCreateWindowCallback(COLSTAT_create_selection_window, weighted_filter->get_column_stat()));
1421    aws->sens_mask(AWM_EXP);
1422    aws->create_button("SELECT_CSP", AWAR_COLUMNSTAT_NAME);
1423    aws->sens_mask(AWM_ALL);
1424
1425    aws->at("alignment");
1426    awt_create_ALI_selection_list(gb_main, aws, AWAR_ALIGNMENT, "*=");
1427
1428    aws->at("tree");
1429    awt_create_TREE_selection_list(gb_main, aws, AWAR_TREE, true);
1430
1431    aws->callback(makeWindowCallback(pars_start_cb, weighted_filter, cmds));
1432    aws->at("go");
1433    aws->create_button("GO", "GO", "G");
1434
1435    return aws;
1436}
1437
1438static void create_parsimony_variables(AW_root *aw_root, AW_default db) {
1439    // kernighan
1440
1441    aw_root->awar_float("genetic/kh/nodes", 1.7, db);
1442    aw_root->awar_int("genetic/kh/maxdepth", 15, db);
1443    aw_root->awar_int("genetic/kh/incdepth", 5, db);
1444
1445    aw_root->awar_int("genetic/kh/static/enable", 1, db);
1446    aw_root->awar_int("genetic/kh/static/depth0", 2, db);
1447    aw_root->awar_int("genetic/kh/static/depth1", 2, db);
1448    aw_root->awar_int("genetic/kh/static/depth2", 2, db);
1449    aw_root->awar_int("genetic/kh/static/depth3", 2, db);
1450    aw_root->awar_int("genetic/kh/static/depth4", 1, db);
1451
1452    aw_root->awar_int("genetic/kh/dynamic/enable", 1,   db);
1453    aw_root->awar_int("genetic/kh/dynamic/start",  100, db);
1454    aw_root->awar_int("genetic/kh/dynamic/maxx",   6,   db);
1455    aw_root->awar_int("genetic/kh/dynamic/maxy",   150, db);
1456
1457    aw_root->awar_int("genetic/kh/function_type", AP_QUADRAT_START, db);
1458
1459    awt_create_dtree_awars(aw_root, db);
1460}
1461
1462static void pars_create_all_awars(AW_root *awr, AW_default aw_def, GBDATA *gb_main) {
1463    awr->awar_string(AWAR_SPECIES_NAME, "",     gb_main);
1464    awr->awar_string(AWAR_FOOTER,       "",     aw_def);
1465    awr->awar_string(AWAR_FILTER_NAME,  "none", gb_main);
1466
1467    {
1468        GB_transaction  ta(gb_main);
1469        char           *dali = GBT_get_default_alignment(gb_main);
1470
1471        awr->awar_string(AWAR_ALIGNMENT, dali, gb_main)->write_string(dali);
1472        free(dali);
1473    }
1474
1475    awt_set_awar_to_valid_filter_good_for_tree_methods(gb_main, awr, AWAR_FILTER_NAME);
1476
1477    awr->awar_string(AWAR_FILTER_FILTER,    "", gb_main);
1478    awr->awar_string(AWAR_FILTER_ALIGNMENT, "", aw_def);
1479
1480    awr->awar(AWAR_FILTER_ALIGNMENT)->map(AWAR_ALIGNMENT);
1481
1482    awr->awar_int(AWAR_PARS_TYPE, PARS_WAGNER, gb_main);
1483
1484    {
1485        GB_transaction  ta(gb_main);
1486        GBDATA *gb_tree_name = GB_search(gb_main, AWAR_TREE, GB_STRING);
1487        char   *tree_name    = GB_read_string(gb_tree_name);
1488
1489        awr->awar_string(AWAR_TREE, "", aw_def)->write_string(tree_name);
1490        free(tree_name);
1491    }
1492
1493    awr->awar_int(AWAR_PARSIMONY,      0, aw_def);
1494    awr->awar_int(AWAR_BEST_PARSIMONY, 0, aw_def);
1495    awr->awar_int(AWAR_STACKPOINTER,   0, aw_def);
1496
1497    create_parsimony_variables(awr, gb_main);
1498    create_nds_vars(awr, aw_def, gb_main);
1499
1500#if defined(DEBUG)
1501    AWT_create_db_browser_awars(awr, aw_def);
1502#endif // DEBUG
1503
1504    GB_ERROR error = ARB_init_global_awars(awr, aw_def, gb_main);
1505    if (error) aw_message(error);
1506}
1507
1508static AW_root *AD_map_viewer_aw_root = 0;
1509
1510void PARS_map_viewer(GBDATA *gb_species, AD_MAP_VIEWER_TYPE vtype) {
1511    if (vtype == ADMVT_SELECT && AD_map_viewer_aw_root && gb_species) {
1512        AD_map_viewer_aw_root->awar(AWAR_SPECIES_NAME)->write_string(GBT_read_name(gb_species));
1513    }
1514}
1515
1516int ARB_main(int argc, char *argv[]) {
1517    aw_initstatus();
1518
1519    GB_shell shell;
1520    AW_root *aw_root      = AWT_create_root("pars.arb", "ARB_PARS", need_macro_ability(), &argc, &argv);
1521    AD_map_viewer_aw_root = aw_root;
1522
1523    ap_main     = new AP_main;
1524    GLOBAL_PARS = new ArbParsimony();
1525
1526    const char *db_server = ":";
1527
1528    PARS_commands cmds;
1529
1530    while (argc>=2 && argv[1][0] == '-') {
1531        argc--;
1532        argv++;
1533        if (!strcmp(argv[0], "-quit"))                   cmds.quit = 1;
1534        else if (!strcmp(argv[0], "-add_marked"))        cmds.add_marked = 1;
1535        else if (!strcmp(argv[0], "-add_selected"))      cmds.add_selected = 1;
1536        else if (!strcmp(argv[0], "-calc_branchlengths")) cmds.calc_branch_lengths = 1;
1537        else if (!strcmp(argv[0], "-calc_bootstrap"))    cmds.calc_bootstrap = 1;
1538        else {
1539            fprintf(stderr, "Unknown option '%s'\n", argv[0]);
1540
1541            printf("    Options:                Meaning:\n"
1542                   "\n"
1543                   "    -add_marked             add marked species   (without changing topology)\n"
1544                   "    -add_selected           add selected species (without changing topology)\n"
1545                   "    -calc_branchlengths     calculate branch lengths only\n"
1546                   "    -calc_bootstrap         estimate bootstrap values\n"
1547                   "    -quit                   quit after performing operations\n"
1548                   );
1549
1550            exit(EXIT_FAILURE);
1551        }
1552    }
1553
1554
1555    if (argc==2) db_server = argv[1];
1556
1557    GB_ERROR error = ap_main->open(db_server);
1558    if (!error) {
1559        GBDATA *gb_main = ap_main->get_gb_main();
1560        error           = configure_macro_recording(aw_root, "ARB_PARS", gb_main);
1561
1562        if (!error) {
1563#if defined(DEBUG)
1564            AWT_announce_db_to_browser(gb_main, GBS_global_string("ARB-database (%s)", db_server));
1565#endif // DEBUG
1566
1567            pars_create_all_awars(aw_root, AW_ROOT_DEFAULT, gb_main);
1568
1569            AW_window *aww = create_pars_init_window(aw_root, &cmds);
1570            aww->show();
1571
1572            AWT_install_cb_guards();
1573            aw_root->main_loop();
1574        }
1575    }
1576
1577    if (error) aw_popup_exit(error);
1578    return EXIT_SUCCESS;
1579}
1580
1581
1582// --------------------------------------------------------------------------------
1583
1584#ifdef UNIT_TESTS
1585#include <arb_file.h>
1586#include <arb_diff.h>
1587#include <test_unit.h>
1588#include <AP_seq_dna.hxx>
1589#include <AP_seq_protein.hxx>
1590#include "test_env.h"
1591
1592#if defined(DEVEL_RALF)
1593// #define AUTO_UPDATE_IF_CHANGED // uncomment to auto update expected results
1594#endif
1595
1596arb_test::match_expectation topologyEquals(AP_tree_nlen *root_node, const char *treefile_base) {
1597    using namespace   arb_test;
1598    expectation_group fulfilled;
1599
1600    char *outfile  = GBS_global_string_copy("pars/%s.tree", treefile_base);
1601    char *expected = GBS_global_string_copy("%s.expected", outfile);
1602    bool  update   = false;
1603
1604    {
1605        FILE *out    = fopen(outfile, "wt");
1606        fulfilled.add(that(out).does_differ_from_NULL());
1607        char *newick = GBT_tree_2_newick(root_node, NewickFormat(nLENGTH|nWRAP), false);
1608        fputs(newick, out);
1609        free(newick);
1610        fclose(out);
1611    }
1612
1613    if (GB_is_regularfile(expected)) {
1614        bool match_exp_topo = textfiles_have_difflines(outfile,expected,0);
1615#if defined(AUTO_UPDATE_IF_CHANGED)
1616        if (!match_exp_topo) update = true;
1617#endif
1618        if (!update) fulfilled.add(that(match_exp_topo).is_equal_to(true));
1619    }
1620    else {
1621        update = true;
1622    }
1623
1624    if (update) TEST_COPY_FILE(outfile, expected);
1625    TEST_EXPECT_ZERO_OR_SHOW_ERRNO(GB_unlink(outfile));
1626
1627    free(expected);
1628    free(outfile);
1629
1630    return all().ofgroup(fulfilled);
1631}
1632#define TEST_EXPECT_SAVED_TOPOLOGY(env,exp_topo) TEST_EXPECTATION(topologyEquals(env.root_node(), exp_topo))
1633
1634#define TEST_EXPECT_PARSVAL(env,exp_pars)  TEST_EXPECT_SIMILAR(env.root_node()->costs(), exp_pars, 0.001);
1635
1636enum TopoMod {
1637    MOD_REMOVE_MARKED,
1638
1639    MOD_QUICK_ADD,
1640    MOD_ADD_NNI,
1641
1642    MOD_ADD_PARTIAL,
1643
1644    MOD_CALC_LENS,
1645    MOD_OPTI_NNI,
1646    MOD_OPTI_GLOBAL,
1647};
1648
1649template <typename SEQ>
1650static void modifyTopology(PARSIMONY_testenv<SEQ>& env, TopoMod mod) {
1651    switch (mod) {
1652        case MOD_REMOVE_MARKED:
1653            env.graphic_tree()->get_tree_root()->remove_leafs(AWT_RemoveType(AWT_REMOVE_BUT_DONT_FREE|AWT_REMOVE_MARKED));
1654            break;
1655
1656        case MOD_QUICK_ADD:
1657            nt_reAdd(env.graphic_tree(), NT_ADD_MARKED, true);
1658            break;
1659
1660        case MOD_ADD_NNI:
1661            nt_reAdd(env.graphic_tree(), NT_ADD_MARKED, false);
1662            break;
1663
1664        case MOD_ADD_PARTIAL:
1665            nt_add_partial(env.graphic_tree());
1666            break;
1667
1668        case MOD_CALC_LENS:
1669            calc_branchlengths(env.graphic_tree());
1670            break;
1671
1672        case MOD_OPTI_NNI:
1673            recursiveNNI(env.graphic_tree());
1674            break;
1675
1676        case MOD_OPTI_GLOBAL:
1677            optimizeTree(env.graphic_tree());
1678            break;
1679    }
1680}
1681
1682template <typename SEQ>
1683static arb_test::match_expectation modifyingTopoResultsIn(TopoMod mod, const char *topo, int pars_expected, PARSIMONY_testenv<SEQ>& env, bool restore) {
1684    using namespace   arb_test;
1685    expectation_group fulfilled;
1686
1687    TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
1688
1689    Level upc = env.get_user_push_counter();
1690    Level fl  = env.get_frame_level();
1691
1692    if (restore) {
1693        env.push();
1694        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
1695    }
1696
1697    modifyTopology(env, mod);
1698    if (topo) fulfilled.add(topologyEquals(env.root_node(), topo));
1699    if (pars_expected != -1) {
1700        fulfilled.add(that(env.root_node()->costs()).fulfills(epsilon_similar(0.001), pars_expected));
1701    }
1702
1703    if (restore) {
1704        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
1705        env.pop();
1706    }
1707
1708    TEST_EXPECT_EQUAL(fl, env.get_frame_level());
1709    TEST_EXPECT_EQUAL(upc, env.get_user_push_counter());
1710
1711    TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
1712
1713    return all().ofgroup(fulfilled);
1714}
1715
1716static GBDATA *copy_to(GBDATA *gb_species, const char *newShortname) {
1717    GBDATA *gb_species_data = GB_get_father(gb_species);
1718    GBDATA *gb_new_species  = GB_create_container(gb_species_data, "species");
1719
1720    GB_ERROR error = NULL;
1721    if (!gb_new_species) {
1722        error = GB_await_error();
1723    }
1724    else {
1725        error = GB_copy(gb_new_species, gb_species);
1726        if (!error) error = GBT_write_string(gb_new_species, "name", newShortname);
1727    }
1728
1729    ap_assert(contradicted(gb_new_species, error));
1730    return gb_new_species;
1731}
1732
1733inline void mark_only(GBDATA *gb_species) {
1734    GBDATA         *gb_main = GB_get_root(gb_species);
1735    GB_transaction  ta(gb_main);
1736    GBT_mark_all(gb_main, 0);
1737    GB_write_flag(gb_species, 1);
1738}
1739inline void mark(GBDATA *gb_species) {
1740    GBDATA         *gb_main = GB_get_root(gb_species);
1741    GB_transaction  ta(gb_main);
1742    GB_write_flag(gb_species, 1);
1743}
1744inline void mark_all(GBDATA *gb_main) {
1745    GB_transaction  ta(gb_main);
1746    GBT_mark_all(gb_main, 0);
1747}
1748
1749inline int is_partial(GBDATA *gb_species) {
1750    GB_transaction ta(gb_species);
1751    return GBT_is_partial(gb_species, -1, false);
1752}
1753
1754template <typename SEQ>
1755static arb_test::match_expectation addedAsBrotherOf(const char *name, const char *allowedBrothers, PARSIMONY_testenv<SEQ>& env) {
1756    using namespace   arb_test;
1757    expectation_group fulfilled;
1758
1759    AP_tree_nlen *node_in_tree = env.root_node()->findLeafNamed(name);
1760    fulfilled.add(that(node_in_tree).does_differ_from(NULL));
1761
1762    const char *brother = node_in_tree->get_brother()->name;
1763    fulfilled.add(that(allowedBrothers).does_contain(brother));
1764
1765    return all().ofgroup(fulfilled);
1766}
1767
1768template <typename SEQ>
1769static arb_test::match_expectation addingPartialResultsIn(GBDATA *gb_added_species, const char *allowedBrothers, const char *topo, int pars_expected, PARSIMONY_testenv<SEQ>& env) {
1770    using namespace   arb_test;
1771    expectation_group fulfilled;
1772
1773    mark_only(gb_added_species);
1774    fulfilled.add(modifyingTopoResultsIn(MOD_ADD_PARTIAL, topo, pars_expected, env, false));
1775    fulfilled.add(that(is_partial(gb_added_species)).is_equal_to(1));
1776
1777    const char *name = GBT_read_name(gb_added_species);
1778    fulfilled.add(addedAsBrotherOf(name, allowedBrothers, env));
1779
1780    return all().ofgroup(fulfilled);
1781}
1782
1783static GBDATA *createPartialSeqFrom(GBDATA *gb_main, const char *aliname, const char *dest_species, const char *source_species, int startPos, int endPos) {
1784    GB_transaction ta(gb_main);
1785
1786    GBDATA *gb_result         = NULL;
1787    GBDATA *gb_source_species = GBT_expect_species(gb_main, source_species);
1788
1789    if (gb_source_species) {
1790        GBDATA *gb_dest_species = copy_to(gb_source_species, dest_species);
1791        GBDATA *gb_dest_seq     = GBT_find_sequence(gb_dest_species, aliname); // =same as source seq
1792        char   *seq             = GB_read_string(gb_dest_seq);
1793
1794        if (seq) {
1795            int maxPos = strlen(seq)-1;
1796
1797            startPos = std::min(startPos, maxPos);
1798            endPos   = std::min(endPos, maxPos);
1799
1800            if (startPos>0) memset(seq, '.', startPos);
1801            if (endPos<maxPos) memset(seq+endPos+1, '.', maxPos-endPos);
1802
1803            GB_ERROR error     = GB_write_string(gb_dest_seq, seq);
1804            if (error) GB_export_error(error);
1805            else {
1806                gb_result = gb_dest_species; // success
1807#if defined(DEBUG)
1808                fprintf(stderr, "created partial '%s' from '%s' (seq='%s')\n", dest_species, source_species, seq);
1809#endif
1810            }
1811
1812            free(seq);
1813        }
1814    }
1815
1816    return gb_result;
1817}
1818
1819static GB_ERROR modifyOneBase(GBDATA *gb_species, const char *aliname, char cOld, char cNew) {
1820    GB_transaction ta(gb_species);
1821    GB_ERROR       error = "failed to modifyOneBase";
1822
1823    GBDATA *gb_seq = GBT_find_sequence(gb_species, aliname);
1824    if (gb_seq) {
1825        char *seq = GB_read_string(gb_seq);
1826        if (seq) {
1827            char *B = strchr(seq, cOld);
1828            if (!B) {
1829                error = "does not contain base in modifyOneBase";
1830            }
1831            else {
1832                B[0]  = cNew;
1833                error = GB_write_string(gb_seq, seq);
1834            }
1835            free(seq);
1836        }
1837    }
1838
1839    return error;
1840}
1841
1842void TEST_nucl_tree_modifications() {
1843    const char *aliname = "ali_5s";
1844
1845    PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", aliname);
1846    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
1847    TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-initial");
1848
1849    const int PARSIMONY_ORG = 301;
1850    TEST_EXPECT_EQUAL(env.combines_performed(), 0);
1851    TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
1852    TEST_EXPECT_EQUAL(env.combines_performed(), 14);
1853
1854    // Note: following code leaks father nodes and edges
1855    // suppressed in valgrind via ../SOURCE_TOOLS/arb.supp@TEST_nucl_tree_modifications
1856
1857    // [NUCOPTI] opposed to protein tests below the initial tree here is NOT optimized! compare .@PROTOPTI
1858    // -> removing and adding species produces a better tree
1859    //
1860    // diff initial->removed  : http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/nucl-removed.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fnucl-initial.tree.expected
1861    // diff initial->add-quick: http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/nucl-add-quick.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fnucl-initial.tree.expected
1862    // diff initial->add-NNI:   http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/nucl-add-NNI.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fnucl-initial.tree.expected
1863
1864    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, "nucl-removed",   PARSIMONY_ORG-93, env, true)); // test remove-marked only (same code as part of nt_reAdd)
1865    TEST_EXPECT_EQUAL(env.combines_performed(), 3);
1866
1867    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD,     "nucl-add-quick", PARSIMONY_ORG-23, env, true)); // test quick-add
1868    TEST_EXPECT_EQUAL(env.combines_performed(), 591);
1869
1870    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_ADD_NNI,       "nucl-add-NNI",   PARSIMONY_ORG-25, env, true)); // test add + NNI
1871    TEST_EXPECT_EQUAL(env.combines_performed(), 925);
1872
1873    // @@@ test optimize etc.
1874
1875    // test partial-add
1876    {
1877        GBDATA *gb_main = env.gbmain();
1878
1879        // create 2 non-overlapping partial species
1880        const int  SPLIT    = 55;
1881        GBDATA    *CorGlutP = createPartialSeqFrom(gb_main, aliname, "CorGlutP", "CorGluta", 0,       SPLIT);
1882        GBDATA    *CloButyP = createPartialSeqFrom(gb_main, aliname, "CloButyP", "CloButyr", SPLIT+1, INT_MAX);
1883        GBDATA    *CloButyM = createPartialSeqFrom(gb_main, aliname, "CloButyM", "CloButyr", SPLIT+1, INT_MAX);
1884        TEST_EXPECT_NO_ERROR(modifyOneBase(CloButyM, aliname, 'G', 'C')); // change first 'G' into 'C'
1885
1886        // add CloButyP (and undo)
1887        {
1888            env.push();
1889            // CloButyr and CloButy2 do not differ in seq-range of partial -> any of both may be chosen as brother.
1890            // behavior should be changed with #605
1891            TEST_EXPECTATION(addingPartialResultsIn(CloButyP, "CloButyr;CloButy2", "nucl-addPart-CloButyP", PARSIMONY_ORG, env));
1892            TEST_EXPECT_EQUAL(env.combines_performed(), 6);
1893            env.pop();
1894        }
1895
1896        {
1897            env.push();
1898            TEST_EXPECTATION(addingPartialResultsIn(CorGlutP, "CorGluta",          "nucl-addPart-CorGlutP",          PARSIMONY_ORG, env)); // add CorGlutP
1899            TEST_EXPECT_EQUAL(env.combines_performed(), 5); // @@@ partial-add should not perform combines at all
1900            TEST_EXPECTATION(addingPartialResultsIn(CloButyP, "CloButyr;CloButy2", "nucl-addPart-CorGlutP-CloButyP", PARSIMONY_ORG, env)); // also add CloButyP
1901            TEST_EXPECT_EQUAL(env.combines_performed(), 6);
1902            env.pop();
1903        }
1904
1905        // now add CorGlutP as full, then CloButyP and CloButyM as partials
1906        {
1907            env.push();
1908
1909            mark_only(CorGlutP);
1910            {
1911                GB_transaction  ta(gb_main);
1912                TEST_EXPECT_NO_ERROR(GBT_write_int(CorGlutP, "ARB_partial", 0)); // revert species to "full"
1913            }
1914
1915            TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD, "nucl-addPartialAsFull-CorGlutP", PARSIMONY_ORG, env, false));
1916            TEST_EXPECT_EQUAL(env.combines_performed(), 254);
1917            TEST_EXPECT_EQUAL(is_partial(CorGlutP), 0); // check CorGlutP was added as full sequence
1918            TEST_EXPECTATION(addedAsBrotherOf("CorGlutP", "CorGluta", env)); // partial created from CorGluta gets inserted next to CorGluta
1919
1920            // add CloButyP as partial.
1921            // as expected it is placed next to matching full sequences (does not differ in partial range)
1922            TEST_EXPECTATION(addingPartialResultsIn(CloButyP, "CloButyr;CloButy2", NULL, PARSIMONY_ORG, env));
1923            TEST_EXPECT_EQUAL(env.combines_performed(), 6);
1924
1925            // CloButyM differs slightly in overlap with CloButyr/CloButy2, but has no overlap with CorGlutP
1926            // reproduces bug described in #609
1927            TEST_EXPECTATION(addingPartialResultsIn(CloButyM, "CorGlutP", "nucl-addPart-bug609",
1928                                                    PARSIMONY_ORG+9, // @@@ known bug - partial should not affect parsimony value; possibly related to ../HELP_SOURCE/oldhelp/pa_partial.hlp@WARNINGS
1929                                                    env));
1930            TEST_EXPECT_EQUAL(env.combines_performed(), 6);
1931            env.pop();
1932        }
1933    }
1934
1935#if 0
1936    // @@@ crashes here with missing sequence (caused by pop()), but works below in TEST_optimizations
1937
1938    // test branchlength calculation
1939    // (optimizations below implicitely recalculate branchlengths)
1940    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_CALC_LENS, "nucl-calclength", PARSIMONY_ORG-17, env, false));
1941#endif
1942}
1943
1944void TEST_optimizations_some() {
1945    const char *aliname = "ali_5s";
1946
1947    PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", aliname);
1948    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
1949    TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-initial");
1950
1951    const int PARSIMONY_ORG = 301;
1952    TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
1953    TEST_EXPECT_EQUAL(env.combines_performed(), 14);
1954
1955    // test branchlength calculation
1956    // (optimizations below implicitely recalculate branchlengths)
1957    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_CALC_LENS, "nucl-calclength", PARSIMONY_ORG, env, false));
1958    TEST_EXPECT_EQUAL(env.combines_performed(), 142);
1959
1960    // test optimize (some)
1961    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_NNI, "nucl-opti-NNI", PARSIMONY_ORG-17, env, true)); // test recursive NNI
1962    TEST_EXPECT_EQUAL(env.combines_performed(), 581);
1963
1964    const int      PARSIMONY_OPTI = 272; // may depend on seed
1965    const unsigned seed           = 1417001558;
1966    GB_random_seed(seed);
1967    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "nucl-opti-global", PARSIMONY_OPTI, env, true)); // test recursive NNI+KL
1968    TEST_EXPECT_EQUAL(env.combines_performed(), 40074);
1969
1970#if 0
1971    // @@@ test results changed by adding the MOD_OPTI_GLOBAL test above (rel #620)
1972    // test optimize (all)
1973    mark_all(env.gbmain());
1974    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_NNI, "nucl-opti-NNI", PARSIMONY_ORG-17, env, true)); // test recursive NNI
1975
1976    // in a fresh environment it works as before (see TEST_optimizations_all)
1977#endif
1978}
1979
1980void TEST_optimizations_all() {
1981    const char *aliname = "ali_5s";
1982
1983    PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", aliname);
1984    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
1985    TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-initial");
1986
1987    const int PARSIMONY_ORG = 301;
1988    TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
1989    TEST_EXPECT_EQUAL(env.combines_performed(), 14);
1990
1991    // test optimize (all)
1992    mark_all(env.gbmain());
1993    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_NNI, "nucl-opti-NNI", PARSIMONY_ORG-17, env, true)); // test recursive NNI
1994    TEST_EXPECT_EQUAL(env.combines_performed(), 581);
1995
1996    const int      PARSIMONY_OPTI = 272; // may depend on seed
1997    const unsigned seed           = 1417001558;
1998    GB_random_seed(seed);
1999    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "nucl-opti-global", PARSIMONY_OPTI, env, true)); // test recursive NNI+KL
2000    TEST_EXPECT_EQUAL(env.combines_performed(), 40074); // @@@ there is no difference in number of combines between TEST_optimizations_some and TEST_optimizations_all. why not?
2001}
2002
2003void TEST_prot_tree_modifications() {
2004    const char *aliname = "ali_tuf_pro";
2005
2006    PARSIMONY_testenv<AP_sequence_protein> env("TEST_prot.arb", aliname);
2007    TEST_EXPECT_NO_ERROR(env.load_tree("tree_prot_opti"));
2008    TEST_EXPECT_SAVED_TOPOLOGY(env, "prot-initial");
2009
2010    const int PARSIMONY_ORG = 917;
2011    TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
2012    TEST_EXPECT_EQUAL(env.combines_performed(), 10);
2013
2014    // Note: following code leaks father nodes and edges
2015    // suppressed in valgrind via ../SOURCE_TOOLS/arb.supp@TEST_prot_tree_modifications
2016
2017    // [PROTOPTI] opposed to nucleid tests above the initial tree here is already optimized! compare .@NUCOPTI
2018    // -> adding species approximately reproduces initial topology
2019    //
2020    // diff initial->add-quick: http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/prot-add-quick.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fprot-initial.tree.expected
2021    // diff initial->add-NNI:   http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/prot-add-NNI.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fprot-initial.tree.expected
2022    //
2023    // Note: comparing these two diffs also demonstrates why quick-adding w/o NNI suffers
2024
2025    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, "prot-removed",   PARSIMONY_ORG-123, env, true)); // test remove-marked only (same code as part of nt_reAdd)
2026    TEST_EXPECT_EQUAL(env.combines_performed(), 5);
2027    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD,     "prot-add-quick", PARSIMONY_ORG+1,   env, true)); // test quick-add
2028    TEST_EXPECT_EQUAL(env.combines_performed(), 302);
2029    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_ADD_NNI,       "prot-add-NNI",   PARSIMONY_ORG,     env, true)); // test add + NNI
2030    TEST_EXPECT_EQUAL(env.combines_performed(), 471);
2031
2032    // @@@ test optimize etc.
2033
2034    // test partial-add
2035    {
2036        GBDATA *gb_main = env.gbmain();
2037
2038        // create 2 non-overlapping partial species
2039        GBDATA    *MucRaceP = createPartialSeqFrom(gb_main, aliname, "MucRaceP", "MucRacem", 0, 60);
2040        GBDATA    *StrCoelP = createPartialSeqFrom(gb_main, aliname, "StrCoelP", "StrCoel9", 66-1, 184-1);
2041        GBDATA    *StrCoelM = createPartialSeqFrom(gb_main, aliname, "StrCoelM", "StrCoel9", 66-1, 184-1);
2042        TEST_EXPECT_NO_ERROR(modifyOneBase(StrCoelM, aliname, 'Y', 'H')); // change first 'Y' into 'H'
2043
2044        // add StrCoelP (and undo)
2045        {
2046            env.push();
2047            // StrCoel9 and StrRamo3 do not differ in seq-range of partial -> any of both may be chosen as brother.
2048            // behavior should be changed with #605
2049            // TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "StrCoel9;StrRamo3", "prot-addPart-StrCoelP", PARSIMONY_ORG+114, env));
2050            TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "AbdGlauc", "prot-addPart-StrCoelP", PARSIMONY_ORG+114, env)); // @@@ inserted completely wrong?
2051            TEST_EXPECT_EQUAL(env.combines_performed(), 4);
2052            env.pop();
2053        }
2054
2055        {
2056            env.push();
2057            TEST_EXPECTATION(addingPartialResultsIn(MucRaceP, "MucRacem",          "prot-addPart-MucRaceP",          PARSIMONY_ORG+7, env)); // add MucRaceP
2058            TEST_EXPECT_EQUAL(env.combines_performed(), 6);
2059            // TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "StrCoel9;StrRamo3", "prot-addPart-MucRaceP-StrCoelP", PARSIMONY_ORG+114, env)); // also add StrCoelP
2060            TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "AbdGlauc", "prot-addPart-MucRaceP-StrCoelP", PARSIMONY_ORG+114+7, env)); // also add StrCoelP // @@@ same misplacement as above
2061            TEST_EXPECT_EQUAL(env.combines_performed(), 4);
2062            env.pop();
2063        }
2064
2065        // now add MucRaceP as full, then StrCoelP and StrCoelM as partials
2066        {
2067            env.push();
2068
2069            mark_only(MucRaceP);
2070            {
2071                GB_transaction  ta(gb_main);
2072                TEST_EXPECT_NO_ERROR(GBT_write_int(MucRaceP, "ARB_partial", 0)); // revert species to "full"
2073            }
2074
2075            TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD, "prot-addPartialAsFull-MucRaceP", PARSIMONY_ORG+7, env, false));
2076            TEST_EXPECT_EQUAL(env.combines_performed(), 177);
2077            TEST_EXPECT_EQUAL(is_partial(MucRaceP), 0); // check CorGlutP was added as full sequence
2078            // TEST_EXPECTATION(addedAsBrotherOf("MucRaceP", "MucRacem", env)); // partial created from MucRacem gets inserted next to MucRacem
2079            TEST_EXPECTATION(addedAsBrotherOf("MucRaceP", "AbdGlauc", env)); // partial created from MucRacem gets inserted next to MucRacem // @@@ misplacement if added as full sequence (as partial placement is ok)
2080
2081            // add StrCoelP as partial.
2082            // as expected it is placed next to matching full sequences (does not differ in partial range)
2083            // TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "StrCoel9;StrRamo3", NULL, PARSIMONY_ORG, env));
2084            TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "MucRaceP", NULL, PARSIMONY_ORG+114+2, env)); // @@@ same position as when adding it partial
2085            TEST_EXPECT_EQUAL(env.combines_performed(), 5);
2086
2087            // StrCoelM differs slightly in overlap with StrCoel9/StrRamo3, but has no overlap with MucRaceP
2088            // reproduces bug described in #609
2089            // @@@ does not demonstrate anything atm (because StrCoelP already is placed next to MucRaceP above)
2090            TEST_EXPECTATION(addingPartialResultsIn(StrCoelM, "StrCoelP", "prot-addPart-bug609",
2091                                                    PARSIMONY_ORG+114+3, // @@@ known bug - partial should not affect parsimony value; possibly related to ../HELP_SOURCE/oldhelp/pa_partial.hlp@WARNINGS
2092                                                    env));
2093            TEST_EXPECT_EQUAL(env.combines_performed(), 6);
2094            env.pop();
2095        }
2096    }
2097}
2098
2099void TEST_node_stack() {
2100    // test was used to fix #620
2101
2102    const char *aliname = "ali_5s";
2103
2104    // Note: following code leaks father nodes and edges
2105    // suppressed in valgrind via ../SOURCE_TOOLS/arb.supp@TEST_node_stack
2106
2107    PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", aliname);
2108    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
2109    TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-initial");
2110
2111    const int PARSIMONY_ORG = 301;
2112    TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
2113    TEST_EXPECT_EQUAL(env.combines_performed(), 14);
2114
2115    TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2116
2117    // test set root to CytAquat + pop (works)
2118    {
2119        env.push();
2120        env.root_node()->findLeafNamed("CytAquat")->set_root();
2121        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2122        env.pop();
2123        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2124    }
2125
2126    // test set root to CloButyr + pop (works)
2127    {
2128        env.push();
2129        env.root_node()->findLeafNamed("CloButyr")->set_root();
2130        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2131        env.pop();
2132        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2133    }
2134
2135    // test set root to CloBifer + set root to CloTyrob + pop (works)
2136    // Note: both species are in same subtree (of root)
2137    {
2138        env.push();
2139
2140        env.root_node()->findLeafNamed("CloBifer")->set_root();
2141        env.root_node()->findLeafNamed("CloTyrob")->set_root();
2142
2143        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2144        env.pop();
2145        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2146    }
2147
2148    // test set root to CytAquat + set root to CloButyr + pop (failed, fixed by [13138])
2149    TEST_EXPECT_EQUAL(env.combines_performed(), 0);
2150    for (int calcCostsBetween = 0; calcCostsBetween<2; ++calcCostsBetween) {
2151        TEST_ANNOTATE(GBS_global_string("calcCostsBetween=%i", calcCostsBetween));
2152
2153        TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
2154
2155        env.push();
2156
2157        env.root_node()->findLeafNamed("CytAquat")->set_root();
2158
2159        if (calcCostsBetween) {
2160            TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
2161            TEST_EXPECT_EQUAL(env.combines_performed(), 2);
2162        }
2163
2164        env.root_node()->findLeafNamed("CloButyr")->set_root();
2165
2166        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2167        TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
2168        TEST_EXPECT_EQUAL(env.combines_performed(), 6);
2169
2170        env.pop();
2171
2172        TEST_EXPECT(env.graphic_tree()->get_root_node()->sequence_state_valid());
2173        TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
2174        TEST_EXPECT_EQUAL(env.combines_performed(), 0);
2175        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2176    }
2177    {
2178        env.push();
2179        {
2180            env.push();
2181
2182            env.root_node()->findLeafNamed("CloInnoc")->moveNextTo(env.root_node()->findLeafNamed("CytAquat"), 0.5);
2183            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2184            env.root_node()->findLeafNamed("CloInnoc")->set_root();
2185            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2186            env.root_node()->findLeafNamed("CytAquat")->moveNextTo(env.root_node()->findLeafNamed("CloPaste"), 0.5);
2187            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2188            env.root_node()->findLeafNamed("CloPaste")->set_root();
2189            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2190            env.root_node()->findLeafNamed("CloPaste")->moveNextTo(env.root_node()->findLeafNamed("CloInnoc"), 0.5);
2191            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2192            {
2193                AP_tree_nlen *son_of_brother;
2194                AP_tree_nlen *brother_of_father;
2195
2196                // COVER1: son of root -> grandson of root
2197                {
2198                    AP_tree_nlen *son_of_root = env.graphic_tree()->get_root_node()->get_leftson();
2199                    ap_assert(son_of_root);
2200
2201                    son_of_brother = son_of_root->get_brother()->get_leftson();
2202                    son_of_root->moveNextTo(son_of_brother, 0.5);
2203                    TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2204                }
2205
2206                // COVER2: grandson of root -> son of brother
2207                {
2208                    AP_tree_nlen *son_of_root      = env.graphic_tree()->get_root_node()->get_leftson();
2209                    AP_tree_nlen *grandson_of_root = son_of_root->get_brother()->get_rightson();
2210                    ap_assert(grandson_of_root);
2211
2212                    son_of_brother = grandson_of_root->get_brother()->get_leftson();
2213                    grandson_of_root->moveNextTo(son_of_brother, 0.5);
2214                    TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2215                }
2216
2217                AP_tree_nlen *some_leaf = env.root_node()->findLeafNamed("CloBifer");
2218                ap_assert(some_leaf);
2219
2220                // COVER3: some leaf -> son of brother
2221                son_of_brother = some_leaf->get_brother()->get_leftson();
2222                some_leaf->moveNextTo(son_of_brother, 0.5);
2223                TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2224
2225                // COVER4: some leaf -> son of brother
2226                brother_of_father = some_leaf->get_father()->get_brother();
2227                some_leaf->moveNextTo(brother_of_father, 0.5);
2228                TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2229
2230                // COVER5: move to father
2231                some_leaf->moveNextTo(some_leaf->get_father(), 0.5);
2232                TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2233
2234#if 0
2235                // fails assert in AP_tree::moveNextTo ("already there"). ok!
2236                some_leaf->moveNextTo(some_leaf->get_brother(), 0.5);
2237                TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2238#endif
2239            }
2240
2241            env.pop();
2242            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2243        }
2244        env.pop();
2245        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2246    }
2247
2248    // remove + quick add marked + pop() both works
2249    TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD,     "nucl-add-quick", PARSIMONY_ORG-23, env, true)); // test quick-add
2250
2251    // remove + quick-add marked + pop() quick-add -> corrupts tree
2252    // (root-edge is lost)
2253    {
2254        env.push();
2255        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2256        TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, NULL, -1, env, false)); // test remove-marked only (same code as part of nt_reAdd)
2257        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2258
2259        env.push();
2260        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2261        TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD, NULL, -1, env, false)); // test quick-add (same code as part of nt_reAdd)
2262        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2263        env.pop();
2264        TEST_EXPECT__BROKEN(env.graphic_tree()->get_root_node()->has_valid_edges()); // @@@ doing pop() after quick-adding produces an invalid tree (root-edge missing)
2265
2266        env.pop();
2267        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2268    }
2269
2270    // same as above, but with only 1 species marked
2271    struct {
2272        const char *name;
2273        bool        correct[2];
2274        bool        correct_after_2nd_pop[2];
2275
2276    } testSingle[] = {
2277        { "CytAquat", { false, true  }, { true,  true  } },  // CytAquat is the only grandson of root (CytAquat located in lower subtree)
2278        { "CloBifer", { false, true  }, { false, true  } },  // two father nodes between CloBifer and root (CloBifer located in upper subtree)
2279        { "CloPaste", { false, true  }, { false, true  } },  // two father nodes between CloPaste and root (CloPaste located in upper subtree)
2280        { "CorGluta", { false,  true }, { false, true  } },  // three father nodes between CorGluta and root (CorGluta located in lower subtree); @@@ might be a different problem
2281        { "CelBiazo", { true,  false }, { true,  true  } },  // two father nodes between CelBiazo and root; @@@ might be a different problem
2282
2283        { NULL, { true, true }, { true, true } }
2284    };
2285
2286    for (int i = 0; testSingle[i].name; ++i) {
2287        for (int swapped = 0; swapped<2; ++swapped) {
2288            TEST_ANNOTATE(GBS_global_string("single=%s swapped=%i", testSingle[i].name, swapped));
2289
2290            env.push();
2291            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2292            {
2293                AP_tree_nlen *old_rightson = env.root_node()->get_rightson();
2294                env.root_node()->get_leftson()->get_rightson()->set_root();
2295                old_rightson->get_leftson()->set_root();
2296                old_rightson->set_root();
2297
2298                ap_assert(env.root_node()->get_rightson() == old_rightson);
2299            }
2300            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2301
2302            mark_only(env.root_node()->findLeafNamed(testSingle[i].name)->gb_node);
2303
2304            env.push();
2305            if (swapped) {
2306                env.root_node()->swap_sons();
2307            }
2308
2309            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2310            TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, NULL, -1, env, false)); // test remove-marked only (same code as part of nt_reAdd)
2311            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2312
2313            env.push();
2314            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2315            TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD, NULL, -1, env, false)); // test quick-add (same code as part of nt_reAdd)
2316            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2317            env.pop();
2318
2319            if (testSingle[i].correct[swapped]) {
2320                TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2321            }
2322            else {
2323                TEST_EXPECT__BROKEN(env.graphic_tree()->get_root_node()->has_valid_edges()); // @@@ doing pop() after quick-adding produces an invalid tree (root-edge missing)
2324            }
2325
2326            env.pop();
2327
2328            if (testSingle[i].correct_after_2nd_pop[swapped]) {
2329                TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2330            }
2331            else {
2332                TEST_EXPECT__BROKEN(env.graphic_tree()->get_root_node()->has_valid_edges()); // @@@ doing 2nd pop() (=undo remove) produces an invalid tree
2333            }
2334
2335            env.pop();
2336            TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2337        }
2338    }
2339
2340    // similar to above (remove+add a grandson of root; here grandson is a subtree with 4 species)
2341
2342    for (int remove_from_lower_subtree = 0; remove_from_lower_subtree<2; ++remove_from_lower_subtree) {
2343        TEST_ANNOTATE(GBS_global_string("remove_from_lower_subtree=%i", remove_from_lower_subtree));
2344
2345        // mark a complete subtree (which - as a whole - forms a grandson of root). subtree is located in upper part of the tree
2346        mark_only(env.root_node()->findLeafNamed("CloButy2")->gb_node);
2347        mark(env.root_node()->findLeafNamed("CloButyr")->gb_node);
2348        mark(env.root_node()->findLeafNamed("CloCarni")->gb_node);
2349        mark(env.root_node()->findLeafNamed("CloPaste")->gb_node);
2350
2351        env.push();
2352        if (remove_from_lower_subtree) {
2353            env.root_node()->swap_sons();
2354        }
2355        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2356        TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, NULL, -1, env, false)); // test remove-marked only (same code as part of nt_reAdd)
2357        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2358
2359        env.push();
2360        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2361        TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD, NULL, -1, env, false)); // test quick-add (same code as part of nt_reAdd)
2362        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2363        env.pop();
2364
2365        if (remove_from_lower_subtree) {
2366            TEST_EXPECT__BROKEN(env.graphic_tree()->get_root_node()->has_valid_edges()); // @@@ broken if removed from lower part of tree
2367        }
2368        else {
2369            TEST_EXPECT(env.graphic_tree()->get_root_node()->has_valid_edges()); // @@@ NOT broken if removed from upper part of tree
2370        }
2371
2372        env.pop();
2373        TEST_ASSERT_VALID_TREE(env.graphic_tree()->get_root_node());
2374    }
2375
2376    TEST_EXPECT_EQUAL(env.combines_performed(), 4438); // @@@ distribute
2377}
2378
2379#endif // UNIT_TESTS
2380
2381// --------------------------------------------------------------------------------
2382
Note: See TracBrowser for help on using the repository browser.