source: branches/species/PARSIMONY/PARS_dtree.cxx

Last change on this file was 19613, checked in by westram, 2 months ago
  • reintegrates 'lib' into 'trunk'
    • replace dynamic library AWT by several static libraries: APP, ARB_SPEC, MASKS, CANVAS, MAPKEY, GUI_TK
    • now also check wrong library dependencies for untested units (only4me)
  • adds: log:branches/lib@19578:19612
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.9 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : PARS_dtree.cxx                                    //
4//   Purpose   : phylogenetic tree display                         //
5//               (specialized for ARB_PARSIMONY)                   //
6//                                                                 //
7//   Institute of Microbiology (Technical University Munich)       //
8//   http://www.arb-home.de/                                       //
9//                                                                 //
10// =============================================================== //
11
12#include "PerfMeter.h"
13#include "pars_dtree.hxx"
14#include "pars_main.hxx"
15#include "pars_awars.h"
16#include "ap_tree_nlen.hxx"
17#include "ap_main.hxx"
18
19#include <AP_TreeColors.hxx>
20#include <AP_seq_dna.hxx>
21#include <AP_seq_protein.hxx>
22#include <AP_filter.hxx>
23
24#include <TreeCallbacks.hxx>
25
26#include <ColumnStat.hxx>
27#include <sel_boxes.hxx>
28#include <ali_filter.hxx>
29
30#include <gui_aliview.hxx>
31
32#include <aw_preset.hxx>
33#include <aw_awar.hxx>
34#include <aw_msg.hxx>
35#include <arb_progress.h>
36#include <aw_root.hxx>
37#include <aw_question.hxx>
38
39static void AWT_graphic_parsimony_root_changed(void *cd, AP_tree *old, AP_tree *newroot) {
40    AWT_graphic_tree *agt = (AWT_graphic_tree*)cd;
41    UNCOVERED();
42
43    if (old == agt->get_logical_root()) agt->set_logical_root_to(newroot);
44}
45
46static AliView *pars_generate_aliview(WeightedFilter *pars_weighted_filter) {
47    GBDATA *gb_main = pars_weighted_filter->get_gb_main();
48    char *ali_name;
49    {
50        GB_transaction ta(gb_main);
51        ali_name = GBT_read_string(gb_main, AWAR_ALIGNMENT);
52    }
53    GB_ERROR  error   = NULp;
54    AliView  *aliview = pars_weighted_filter->create_aliview(ali_name, error);
55    if (!aliview) aw_popup_exit(error);
56    free(ali_name);
57    return aliview;
58}
59
60void PARS_tree_init(AWT_graphic_parsimony *agt) {
61    ap_assert(agt->get_root_node());
62    ap_assert(agt == ap_main->get_graphic_tree());
63
64    GBDATA         *gb_main = ap_main->get_gb_main();
65    GB_transaction  ta(gb_main);
66
67    const char *use     = ap_main->get_aliname();
68    long        ali_len = GBT_get_alignment_len(gb_main, use);
69    if (ali_len <= 1) {
70        aw_popup_exit("No valid alignment selected! Try again");
71    }
72
73    agt->get_tree_root()->set_root_changed_callback(AWT_graphic_parsimony_root_changed, agt);
74}
75
76QuadraticThreshold::QuadraticThreshold(KL_DYNAMIC_THRESHOLD_TYPE type, double startx, double maxy, double maxx, double maxDepth, Mutations pars_start) {
77    // set a, b,  and c for quadratic equation y = ax^2 + bx + c
78    switch (type) {
79        default:
80            ap_assert(0);
81            // fall-through
82        case AP_QUADRAT_START:
83            c = startx;
84            a = (startx - maxy) / (maxx * maxx);
85            b = -2.0 * a * maxx;
86            break;
87
88        case AP_QUADRAT_MAX: // unused (experimental)
89            a = - maxy / ((maxDepth -  maxx) * (maxDepth - maxx));
90            b = -2.0 * a * maxx;
91            c = maxy  + a * maxx * maxx;
92            break;
93    }
94    c += pars_start;
95}
96
97void ArbParsimony::kernighan_optimize_tree(AP_tree_nlen *at, const KL_Settings& settings, const Mutations *pars_global_start, bool dumpPerf) {
98    AP_tree_nlen    *oldrootleft  = get_root_node()->get_leftson();
99    AP_tree_nlen    *oldrootright = get_root_node()->get_rightson();
100    Mutations        pars_curr    = get_root_node()->costs();
101    const Mutations  pars_org     = pars_curr;
102
103    OptiPerfMeter performance("KL-recursion", pars_curr);
104
105    // setup KL recursion parameters:
106    KL_params KL;
107    {
108        KL.max_rec_depth = settings.maxdepth; ap_assert(KL.max_rec_depth>0);
109        KL.inc_rec_depth = settings.incdepth;
110        KL.thresFunctor  = QuadraticThreshold(settings.Dynamic.type, settings.Dynamic.start, settings.Dynamic.maxy, settings.Dynamic.maxx, KL.max_rec_depth, pars_curr);
111        KL.rec_type      = KL_RECURSION_TYPE((settings.Dynamic.enabled*AP_DYNAMIK)|(settings.Static.enabled*AP_STATIC));
112
113        for (int x = 0; x<CUSTOM_DEPTHS; ++x) {
114            KL.rec_width[x] = settings.Static.depth[x];
115        }
116        KL.stopAtFoldedGroups = settings.whichEdges&SKIP_FOLDED_EDGES;
117    }
118
119    AP_tree_edge *startEdge     = NULp;
120    AP_tree_nlen *skipNode      = NULp;
121    bool          skipStartEdge = true;
122
123    if (!at->father) {
124        startEdge     = rootEdge();
125        skipStartEdge = false;
126    }
127    else if (at->is_son_of_root()) {
128        startEdge = rootEdge();
129        skipNode  = startEdge->otherNode(at);
130    }
131    else {
132        skipNode  = at->get_father();
133        startEdge = at->edgeTo(skipNode);
134    }
135    ap_assert(startEdge);
136
137    EdgeChain    chain(startEdge, EdgeSpec(SKIP_LEAF_EDGES|settings.whichEdges), true, skipNode, skipStartEdge);
138    arb_progress progress(chain.size());
139
140    if (pars_global_start) {
141        progress.subtitle(GBS_global_string("best=%li (gain=%li)", pars_curr, *pars_global_start-pars_curr));
142    }
143    else {
144        progress.subtitle(GBS_global_string("best=%li", pars_curr));
145    }
146
147    if (skipStartEdge) startEdge->set_visited(true); // avoid traversal beyond startEdge
148
149    while (chain && !progress.aborted()) {
150        AP_tree_edge *edge = *chain; ++chain;
151
152        ap_assert(!edge->is_leaf_edge());
153        ap_assert(implicated(KL.stopAtFoldedGroups, !edge->next_to_folded_group()));
154
155        ap_main->remember();
156
157        bool better_tree_found = edge->kl_rec(KL, 0, pars_curr);
158
159        if (better_tree_found) {
160            ap_main->accept();
161            Mutations pars_new = get_root_node()->costs();
162            KL.thresFunctor.change_parsimony_start(pars_new-pars_curr);
163            pars_curr = pars_new;
164            if (pars_global_start) {
165                progress.subtitle(GBS_global_string("best=%li (gain=%li, KL=%li)", pars_curr, *pars_global_start-pars_curr, pars_org-pars_curr));
166            }
167            else {
168                progress.subtitle(GBS_global_string("best=%li (gain=%li)", pars_curr, pars_org-pars_curr));
169            }
170        }
171        else {
172            ap_main->revert();
173        }
174        progress.inc();
175    }
176
177    if (skipStartEdge) startEdge->set_visited(false); // reallow traversal beyond startEdge
178
179    if (dumpPerf) performance.dump(stdout, pars_curr);
180
181    if (oldrootleft->father == oldrootright) {
182        oldrootleft->set_root();
183    }
184    else {
185        oldrootright->set_root();
186    }
187}
188
189
190
191void ArbParsimony::optimize_tree(AP_tree_nlen *at, const KL_Settings& settings, arb_progress& progress) {
192    AP_tree_nlen    *oldrootleft  = get_root_node()->get_leftson();
193    AP_tree_nlen    *oldrootright = get_root_node()->get_rightson();
194    const Mutations  org_pars     = get_root_node()->costs();
195    Mutations        prev_pars    = org_pars;
196
197    OptiPerfMeter overallPerf("global optimization", org_pars);
198
199    progress.subtitle(GBS_global_string("best=%li", org_pars));
200
201    // define available heuristics
202    enum Heuristic {
203        RECURSIVE_NNI,
204        CUSTOM_KL,
205
206        NO_FURTHER_HEURISTIC,
207        HEURISTIC_COUNT = NO_FURTHER_HEURISTIC,
208    } heuristic = RECURSIVE_NNI;
209
210    struct H_Settings {
211        const char        *name;      // name shown in OptiPerfMeter
212        const KL_Settings *kl;        // ==NULp -> NNI; else KL with these settings
213        bool               repeat;    // retry same heuristic when tree improved
214        Heuristic          onImprove; // continue with this heuristic if improved (repeated or not)
215        Heuristic          onFailure; // continue with this heuristic if NOT improved
216    } heuristic_setting[HEURISTIC_COUNT] = {
217        { "recursive NNIs", NULp,      true,  CUSTOM_KL,     CUSTOM_KL },
218        { "KL-optimizer",   &settings, false, RECURSIVE_NNI, NO_FURTHER_HEURISTIC },
219    };
220
221    Mutations      heu_start_pars = prev_pars;
222    OptiPerfMeter *heuPerf        = NULp;
223
224#if defined(ASSERTION_USED)
225    bool at_is_root = at == rootNode();
226#endif
227
228    {
229        arb_suppress_progress quiet;
230
231        while (heuristic != NO_FURTHER_HEURISTIC && !progress.aborted()) {
232            const H_Settings& hset = heuristic_setting[heuristic];
233            if (!heuPerf) {
234                ap_assert(heu_start_pars == prev_pars);
235                heuPerf = new OptiPerfMeter(hset.name, heu_start_pars);
236            }
237
238            Mutations this_pars(-1);
239            if (hset.kl) {
240                kernighan_optimize_tree(at, *hset.kl, &org_pars, false);
241                this_pars = get_root_node()->costs();
242            }
243            else {
244                this_pars = at->nn_interchange_rec(settings.whichEdges, AP_BL_NNI_ONLY);
245            }
246            ap_assert(this_pars>=0); // ensure this_pars was set
247            ap_assert(this_pars<=prev_pars); // otherwise heuristic worsened the tree
248
249            ap_assert(at_is_root == (at == rootNode()));
250
251            bool      dumpOverall   = false;
252            Heuristic nextHeuristic = heuristic;
253            if (this_pars<prev_pars) { // found better tree
254                prev_pars = this_pars;
255                progress.subtitle(GBS_global_string("best=%li (gain=%li)", this_pars, org_pars-this_pars));
256                if (!hset.repeat) {
257                    nextHeuristic = hset.onImprove;
258                    dumpOverall   = heuristic == CUSTOM_KL;
259                }
260            }
261            else { // last step did not find better tree
262                nextHeuristic = this_pars<heu_start_pars ? hset.onImprove : hset.onFailure;
263            }
264
265            if (nextHeuristic != heuristic) {
266                heuristic      = nextHeuristic;
267                heu_start_pars = this_pars;
268
269                heuPerf->dump(stdout, this_pars);
270                delete heuPerf; heuPerf = NULp;
271            }
272            if (dumpOverall) overallPerf.dumpCustom(stdout, this_pars, "overall (so far)");
273        }
274    }
275
276    if (oldrootleft->father == oldrootright) {
277        oldrootleft->set_root();
278    }
279    else {
280        oldrootright->set_root();
281    }
282
283    overallPerf.dump(stdout, prev_pars);
284}
285
286AWT_graphic_parsimony::AWT_graphic_parsimony(ArbParsimony& parsimony_, GBDATA *gb_main_, AD_map_viewer_cb map_viewer_cb_)
287    : AWT_graphic_tree(AW_root::SINGLETON, gb_main_, map_viewer_cb_),
288      parsimony(parsimony_)
289{}
290
291AP_tree_root *AWT_graphic_parsimony::create_tree_root(AliView *aliview, AP_sequence *seq_prototype, bool insert_delete_cbs) {
292    return new AP_pars_root(aliview, seq_prototype, insert_delete_cbs, &groupScale);
293}
294
295
296void ArbParsimony::generate_tree(WeightedFilter *pars_weighted_filter) {
297    AliView     *aliview   = pars_generate_aliview(pars_weighted_filter);
298    AP_sequence *seq_templ = NULp;
299
300    GBDATA *gb_main = aliview->get_gb_main();
301    {
302        GB_transaction ta(gb_main);
303        bool           is_aa = GBT_is_alignment_protein(gb_main, aliview->get_aliname());
304
305        if (is_aa) seq_templ = new AP_sequence_protein(aliview);
306        else seq_templ       = new AP_sequence_parsimony(aliview);
307    }
308
309    AWT_graphic_parsimony *new_tree = new AWT_graphic_parsimony(*this, aliview->get_gb_main(), PARS_map_viewer);
310    new_tree->init(aliview, seq_templ, true, false);
311    set_tree(new_tree);
312}
313
314AW_gc_manager *AWT_graphic_parsimony::init_devices(AW_window *aww, AW_device *device, AWT_canvas* ntw) {
315    AW_init_color_group_defaults("arb_ntree");
316
317    AW_gc_manager *gc_manager =
318        AW_manage_GC(aww,
319                     ntw->get_gc_base_name(),
320                     device, AWT_GC_MAX, AW_GCM_DATA_AREA,
321                     makeGcChangedCallback(TREE_GC_changed_cb, ntw),
322                     "#AAAA55",
323
324                     // Important note :
325                     // Many gc indices are shared between ABR_NTREE and ARB_PARSIMONY
326                     // e.g. the tree drawing routines use same gc's for drawing both trees
327                     // (keep in sync with ../SL/TREEDISP/TreeDisplay.cxx@init_devices)
328
329                     "Cursor$#FFFFFF",
330                     "Branch remarks$#DBE994",
331                     "+-Bootstrap$#DBE994",    "-B.(limited)$white",
332                     "!1", // reserve GC-number used for "IRS group box" in arb_ntree
333                     "Marked$#FFFF00",
334                     "Some marked$#eeee88",
335                     "Not marked$black",
336                     "Zombies etc.$#cc5924",
337
338                     "!14", // reserve 14 GC-numbers which are used for probe colors in ARB_NTREE
339                            // (namely 'None', 'All' and 'P1' up to 'P12')
340
341                     "&color_groups", // use color groups
342
343                     NULp);
344    return gc_manager;
345}
346
347void AWT_graphic_parsimony::show(AW_device *device) {
348    AP_tree_nlen *root_node = parsimony.get_root_node();
349    AW_awar      *awar_pars = get_root()->awar(AWAR_PARSIMONY);
350    AW_awar      *awar_best = get_root()->awar(AWAR_BEST_PARSIMONY);
351
352    long parsval = root_node ? root_node->costs() : 0;
353    awar_pars->write_int(parsval);
354
355    long best_pars = awar_best->read_int();
356    if (parsval < best_pars || 0==best_pars) awar_best->write_int(parsval);
357
358    AWT_graphic_tree::show(device);
359}
360
361void AWT_graphic_parsimony::handle_command(AW_device *device, AWT_graphic_event& event) {
362    ClickedTarget clicked(this, event.best_click());
363    bool          recalc_branchlengths_on_structure_change = true;
364
365    switch (event.cmd()) {
366        case AWT_MODE_NNI:
367        case AWT_MODE_KERNINGHAN:
368        case AWT_MODE_OPTIMIZE: {
369            const char *what  = NULp;
370            const char *where = NULp;
371
372            switch (event.cmd()) {
373                case AWT_MODE_NNI:        what = "Recursive NNI on"; break;
374                case AWT_MODE_KERNINGHAN: what = "K.L. optimize";    break;
375                case AWT_MODE_OPTIMIZE:   what = "Global optimize";  break;
376                default:                  break;
377            }
378
379            AP_tree_nlen *startNode  = NULp;
380            bool          repeatOpti = true;
381
382            if (event.type() == AW_Mouse_Press) {
383                switch (event.button()) {
384                    case AW_BUTTON_LEFT:
385                        repeatOpti = false;
386                        // fall-through
387                    case AW_BUTTON_RIGHT:
388                        startNode = DOWNCAST(AP_tree_nlen*, clicked.node());
389                        where = (startNode == get_root_node()) ? "tree" : "subtree";
390                        break;
391
392                    default:
393                        break;
394                }
395            }
396
397            if (what && where) {
398                arb_progress progress(GBS_global_string("%s %s", what, where));
399
400                Mutations start_pars = get_root_node()->costs();
401                Mutations curr_pars  = start_pars;
402
403                KL_Settings KL(get_root());
404
405                do {
406                    Mutations prev_pars = curr_pars;
407
408                    switch (event.cmd()) {
409                        case AWT_MODE_NNI:
410                            startNode->nn_interchange_rec(KL.whichEdges, AP_BL_NNI_ONLY);
411                            break;
412                        case AWT_MODE_KERNINGHAN:
413                            parsimony.kernighan_optimize_tree(startNode, KL, &start_pars, true);
414                            break;
415                        case AWT_MODE_OPTIMIZE:
416                            parsimony.optimize_tree(startNode, KL, progress);
417                            repeatOpti = false;   // never loop here (optimize_tree already loops until no further improvement)
418                            break;
419                        default:
420                            repeatOpti = false;
421                            break;
422                    }
423
424                    curr_pars  = get_root_node()->costs();
425                    repeatOpti = repeatOpti && curr_pars<prev_pars;
426                } while (repeatOpti);
427
428                exports.request_save();
429                ASSERT_VALID_TREE(get_root_node());
430            }
431            break;
432        }
433
434        default:
435            recalc_branchlengths_on_structure_change = false;
436            FALLTHROUGH; // unlisted modes trigger branchlength calculation internally when needed
437        case AWT_MODE_MOVE:
438            AWT_graphic_tree::handle_command(device, event);
439            break;
440    }
441
442    if (exports.needs_save() && recalc_branchlengths_on_structure_change) {
443        arb_progress progress("Recalculating branch lengths");
444        rootEdge()->calc_branchlengths();
445        reorderTree(BIG_BRANCHES_TO_TOP); // beautify after recalc_branch_lengths
446    }
447}
448
449
450// --------------------------------------------------------------------------------
451
452#ifdef UNIT_TESTS
453#include <arb_diff.h>
454#include <test_unit.h>
455#include "test_env.h"
456
457template<typename SEQTYPE>
458PARSIMONY_testenv<SEQTYPE>::PARSIMONY_testenv(const char *dbname, const char *aliName)
459    : parsimony(),
460      klSettings(NULp)
461{
462    common_init(dbname);
463    GBDATA         *gb_main = ap_main->get_gb_main();
464    GB_transaction  ta(gb_main);
465
466    size_t aliLength = GBT_get_alignment_len(gb_main, aliName);
467    ap_assert(aliLength>0);
468
469    klSettings = new KL_Settings();
470
471    AP_filter filter(aliLength);
472    if (!filter.is_invalid()) {
473        AP_weights weights(&filter);
474        agt->init(new AliView(gb_main, filter, weights, aliName));
475    }
476}
477
478template PARSIMONY_testenv<AP_sequence_protein>::PARSIMONY_testenv(const char *dbname, const char *aliName);   // explicit instanciation (otherwise link error in unittest)
479template PARSIMONY_testenv<AP_sequence_parsimony>::PARSIMONY_testenv(const char *dbname, const char *aliName); // explicit instanciation (as above, but only happens with gcc 4.6.3/NDEBUG)
480
481
482void TEST_basic_tree_modifications() {
483    PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb");
484    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
485
486    {
487        AP_tree_nlen *root = env.root_node();
488
489        // first check initial state:
490        {
491            AP_tree_members& root_info = root->gr;
492
493            TEST_EXPECT_EQUAL(root_info.grouped,  false);
494            TEST_EXPECT_EQUAL(root_info.hidden,   false);
495            TEST_EXPECT_EQUAL(root_info.mark_sum, 6);
496            TEST_EXPECT_EQUAL(root_info.leaf_sum, 15);
497
498            TEST_EXPECT_SIMILAR(root_info.max_tree_depth, 1.624975, 0.000001);
499            TEST_EXPECT_SIMILAR(root_info.min_tree_depth, 0.341681, 0.000001);
500
501            GB_transaction ta(env.gbmain());
502            GBT_mark_all(env.gbmain(), 0); // unmark all species
503            root->compute_tree();
504            TEST_EXPECT_EQUAL(root_info.mark_sum, 0);
505        }
506
507
508#define B1_TOP "(((((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,CloTyrob:0.009):0.274,CloInnoc:0.371):0.057,CloBifer:0.388):0.124"
509#define B1_BOT "(CloBifer:0.388,(CloInnoc:0.371,(CloTyrob:0.009,(CloTyro2:0.017,(CloTyro3:1.046,CloTyro4:0.061):0.026):0.017):0.274):0.057):0.124"
510#define B2_TOP "(((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010,CloPaste:0.179):0.131"
511#define B2_BOT "(CloPaste:0.179,(CloCarni:0.120,(CloButy2:0.009,CloButyr:0.000):0.564):0.010):0.131"
512
513
514#define B3_LEFT_TOP_SONS "(((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053,CelBiazo:0.059)"
515#define B3_TOP_SONS      B3_LEFT_TOP_SONS ":0.207,CytAquat:0.711"
516#define B3_TOP_SONS_CCR  "((CorAquat:0.187,CorGluta:0.522):0.053,CelBiazo:0.059):0.207,CytAquat:0.711" // CCR = CurCitre removed
517#define B3_TOP           "(" B3_TOP_SONS "):0.081"
518#define B3_BOT           "(CytAquat:0.711,(CelBiazo:0.059,(CorGluta:0.522,(CorAquat:0.084,CurCitre:0.058):0.103):0.053):0.207):0.081"
519
520
521        const char *top_topo    = "((" B1_TOP "," B2_TOP "):0.081," B3_TOP ");";
522        const char *edge_topo   = "((" B1_TOP "," B2_BOT "):0.081," B3_BOT ");";
523        const char *bottom_topo = "(" B3_BOT ",(" B2_BOT "," B1_BOT "):0.081);";
524
525        const char *radial_topo  =
526            "(((CloPaste:0.179,((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010):0.131,"
527            "((CloInnoc:0.371,((CloTyro2:0.017,(CloTyro3:1.046,CloTyro4:0.061):0.026):0.017,CloTyrob:0.009):0.274):0.057,CloBifer:0.388):0.124):0.081,"
528            "((CelBiazo:0.059,((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053):0.207,CytAquat:0.711):0.081);";
529        const char *radial_topo2 =
530            "(((CloBifer:0.388,(CloInnoc:0.371,(((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,CloTyrob:0.009):0.274):0.057):0.124," B2_TOP "):0.081,"
531            "(CytAquat:0.711," B3_LEFT_TOP_SONS ":0.207):0.081);";
532
533        // expect that no mode reproduces another mode:
534        TEST_EXPECT_DIFFERENT(top_topo,    edge_topo);
535        TEST_EXPECT_DIFFERENT(top_topo,    bottom_topo);
536        TEST_EXPECT_DIFFERENT(top_topo,    radial_topo);
537        TEST_EXPECT_DIFFERENT(top_topo,    radial_topo2);
538        TEST_EXPECT_DIFFERENT(edge_topo,   bottom_topo);
539        TEST_EXPECT_DIFFERENT(edge_topo,   radial_topo);
540        TEST_EXPECT_DIFFERENT(edge_topo,   radial_topo2);
541        TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo);
542        TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo2);
543        TEST_EXPECT_DIFFERENT(radial_topo, radial_topo2);
544
545        env.push(); // 1st stack level (=top_topo)
546
547        TEST_EXPECT_VALID_TREE(root);
548
549        TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
550        // test reorder_tree:
551        root->reorder_tree(BIG_BRANCHES_TO_EDGE);     TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo);    env.push(); // 2nd stack level (=edge_topo)
552        root->reorder_tree(BIG_BRANCHES_TO_BOTTOM);   TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo);  env.push(); // 3rd stack level (=bottom_topo)
553        root->reorder_tree(BIG_BRANCHES_TO_CENTER);   TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo);
554        root->reorder_tree(BIG_BRANCHES_ALTERNATING); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo2);
555        root->reorder_tree(BIG_BRANCHES_TO_TOP);      TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
556
557        TEST_EXPECT_VALID_TREE(root);
558
559        // test set root:
560        AP_tree_nlen *CloTyrob = root->findLeafNamed("CloTyrob");
561        TEST_REJECT_NULL(CloTyrob);
562
563        ARB_edge rootEdge(root->get_leftson(), root->get_rightson());
564        CloTyrob->set_root();
565
566        TEST_EXPECT_VALID_TREE(root);
567
568        const char *rootAtCloTyrob_topo =
569            "(CloTyrob:0.004,"
570            "(((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,"
571            "((((" B3_TOP_SONS "):0.162," B2_TOP "):0.124,CloBifer:0.388):0.057,CloInnoc:0.371):0.274):0.004);";
572
573        TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo);
574        env.push(); // 4th stack level (=rootAtCloTyrob_topo)
575
576        TEST_EXPECT_VALID_TREE(root);
577
578        AP_tree_nlen *CelBiazoFather = root->findLeafNamed("CelBiazo")->get_father();
579        TEST_REJECT_NULL(CelBiazoFather);
580        CelBiazoFather->set_root();
581
582        const char *rootAtCelBiazoFather_topo = "(" B3_LEFT_TOP_SONS ":0.104,((" B1_TOP "," B2_TOP "):0.162,CytAquat:0.711):0.104);";
583        TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCelBiazoFather_topo);
584
585        TEST_EXPECT_VALID_TREE(root);
586
587        ARB_edge oldRootEdge(rootEdge.source(), rootEdge.dest());
588        DOWNCAST(AP_tree_nlen*,oldRootEdge.son())->set_root();
589
590        const char *rootSetBack_topo = top_topo;
591        TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo);
592        env.push(); // 5th stack level (=rootSetBack_topo)
593
594        TEST_EXPECT_VALID_TREE(root);
595
596        // test remove:
597        AP_tree_nlen *CurCitre = root->findLeafNamed("CurCitre");
598        TEST_REJECT_NULL(CurCitre);
599        TEST_REJECT_NULL(CurCitre->get_father());
600
601        CurCitre->REMOVE();
602        const char *CurCitre_removed_topo = "((" B1_TOP "," B2_TOP "):0.081,(" B3_TOP_SONS_CCR "):0.081);";
603        // ------------------------------------------------------------------- ^^^ = B3_TOP_SONS minus CurCitre
604        TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo);
605
606        TEST_EXPECT_VALID_TREE(root);
607        TEST_EXPECT_VALID_TREE(CurCitre);
608
609        TEST_EXPECT_EQUAL(root->gr.leaf_sum, 15); // out of date
610        root->compute_tree();
611        TEST_EXPECT_EQUAL(root->gr.leaf_sum, 14);
612
613        env.push(); // 6th stack level (=CurCitre_removed_topo)
614
615        TEST_EXPECT_VALID_TREE(root);
616
617        // test insert:
618        AP_tree_nlen *CloCarni = root->findLeafNamed("CloCarni");
619        TEST_REJECT_NULL(CloCarni);
620        CurCitre->insert(CloCarni); // this creates two extra edges (not destroyed by destroy() below) and one extra node
621
622        const char *CurCitre_inserted_topo = "((" B1_TOP ",(((CloButy2:0.009,CloButyr:0.000):0.564,(CurCitre:0.060,CloCarni:0.060):0.060):0.010,CloPaste:0.179):0.131):0.081,(" B3_TOP_SONS_CCR "):0.081);";
623        TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_inserted_topo);
624
625        TEST_EXPECT_VALID_TREE(root);
626
627        // now check pops:
628        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo); TEST_EXPECT_VALID_TREE(root);
629        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo);      TEST_EXPECT_VALID_TREE(root);
630        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo);   TEST_EXPECT_VALID_TREE(root);
631        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo);           TEST_EXPECT_VALID_TREE(root);
632        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo);             TEST_EXPECT_VALID_TREE(root);
633        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);              TEST_EXPECT_VALID_TREE(root);
634    }
635}
636
637void TEST_calc_bootstraps() {
638    PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", "ali_5s");
639    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
640
641    const char *bs_origi_topo = "(((((((CloTyro3,CloTyro4)'40%',CloTyro2)'0%',CloTyrob)'97%',CloInnoc)'0%',CloBifer)'53%',(((CloButy2,CloButyr),CloCarni)'33%',CloPaste)'97%'),((((CorAquat,CurCitre),CorGluta)'17%',CelBiazo)'40%',CytAquat));";
642    const char *bs_limit_topo = "(((((((CloTyro3,CloTyro4)'87%',CloTyro2)'0%',CloTyrob)'100%',CloInnoc)'87%',CloBifer)'83%',(((CloButy2,CloButyr)'99%',CloCarni)'17%',CloPaste)'56%')'61%',((((CorAquat,CurCitre)'78%',CorGluta)'0%',CelBiazo)'59%',CytAquat)'61%');";
643    const char *bs_estim_topo = "(((((((CloTyro3,CloTyro4)'75%',CloTyro2)'0%',CloTyrob)'100%',CloInnoc)'75%',CloBifer)'78%',(((CloButy2,CloButyr)'99%',CloCarni)'13%',CloPaste)'32%')'53%',((((CorAquat,CurCitre)'74%',CorGluta)'0%',CelBiazo)'56%',CytAquat)'53%');";
644
645    {
646        AP_tree_nlen *root      = env.root_node();
647        AP_tree_edge *root_edge = rootEdge();
648
649        TEST_EXPECT(root && root_edge);
650
651        root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_origi_topo);
652
653        TEST_EXPECT_COMBINES_PERFORMED(env, 0);
654
655        root_edge->nni_rec(ANY_EDGE, AP_BL_MODE(AP_BL_BL_ONLY|AP_BL_BOOTSTRAP_LIMIT),    NULp, true);
656        root->reorder_tree(BIG_BRANCHES_TO_TOP);
657        TEST_EXPECT_NEWICK(nREMARK, root, bs_limit_topo);
658        TEST_EXPECT_COMBINES_PERFORMED(env, 170);
659
660        root_edge->nni_rec(ANY_EDGE, AP_BL_MODE(AP_BL_BL_ONLY|AP_BL_BOOTSTRAP_ESTIMATE), NULp, true);
661        root->reorder_tree(BIG_BRANCHES_TO_TOP);
662        TEST_EXPECT_NEWICK(nREMARK, root, bs_estim_topo);
663        TEST_EXPECT_COMBINES_PERFORMED(env, 156);
664
665        TEST_EXPECT_EQUAL(env.root_node(), root);
666    }
667
668}
669
670void TEST_tree_remove_add_all() {
671    // reproduces crash as described in #527
672    PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", "ali_5s");
673    TEST_EXPECT_NO_ERROR(env.load_tree("tree_nj"));
674
675    const int     LEAFS     = 6;
676    AP_tree_nlen *leaf[LEAFS];
677    const char *name[LEAFS] = {
678        "CloButy2",
679        "CloButyr",
680        "CytAquat",
681        "CorAquat",
682        "CurCitre",
683        "CorGluta",
684    };
685
686    AP_tree_nlen *root = env.root_node();
687
688    for (int i = 0; i<LEAFS; ++i) {
689        leaf[i] = root->findLeafNamed(name[i]);
690        TEST_REJECT_NULL(leaf[i]);
691    }
692
693    TEST_EXPECT_VALID_TREE(root);
694
695    AP_pars_root *troot = leaf[0]->get_tree_root();
696    TEST_REJECT_NULL(troot);
697
698    for (int i = 0; i<LEAFS-1; ++i) {
699        // Note: removing the second to last leaf, "removes" both remaining
700        // leafs (but only destroys their father node)
701
702        TEST_EXPECT_VALID_TREE(root);
703        leaf[i]->REMOVE();
704        TEST_EXPECT_VALID_TREE(leaf[i]);
705    }
706
707    leaf[0]->initial_insert(leaf[1], troot);
708    for (int i = 2; i<LEAFS; ++i) {
709        TEST_EXPECT_VALID_TREE(leaf[i-1]);
710        TEST_EXPECT_VALID_TREE(leaf[i]);
711        leaf[i]->insert(leaf[i-1]);
712    }
713}
714
715#endif // UNIT_TESTS
716
717// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.