source: tags/arb-6.0.3/PARSIMONY/PARS_dtree.cxx

Last change on this file was 12554, checked in by westram, 11 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.5 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : PARS_dtree.cxx                                    //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "pars_dtree.hxx"
12#include "pars_main.hxx"
13#include "pars_debug.hxx"
14#include "ap_tree_nlen.hxx"
15
16#include <AP_seq_dna.hxx>
17#include <AP_seq_protein.hxx>
18#include <AP_filter.hxx>
19
20#include <ColumnStat.hxx>
21#include <awt_sel_boxes.hxx>
22#include <awt_filter.hxx>
23
24#include <gui_aliview.hxx>
25
26#include <aw_preset.hxx>
27#include <aw_awar.hxx>
28#include <aw_msg.hxx>
29#include <arb_progress.h>
30#include <aw_root.hxx>
31#include <aw_question.hxx>
32
33static void AWT_graphic_parsimony_root_changed(void *cd, AP_tree *old, AP_tree *newroot) {
34    AWT_graphic_tree *agt = (AWT_graphic_tree*)cd;
35
36    if (old == agt->displayed_root) agt->displayed_root = newroot;
37}
38
39static AliView *pars_generate_aliview(WeightedFilter *pars_weighted_filter) {
40    GBDATA *gb_main = pars_weighted_filter->get_gb_main();
41    char *ali_name;
42    {
43        GB_transaction ta(gb_main);
44        ali_name = GBT_read_string(gb_main, AWAR_ALIGNMENT);
45    }
46    GB_ERROR  error   = NULL;
47    AliView  *aliview = pars_weighted_filter->create_aliview(ali_name, error);
48    if (!aliview) aw_popup_exit(error);
49    free(ali_name);
50    return aliview;
51}
52
53void PARS_tree_init(AWT_graphic_tree *agt) {
54    ap_assert(agt->get_root_node());
55    ap_assert(agt == ap_main->get_tree_root());
56
57    GB_transaction ta(GLOBAL_gb_main);
58
59    const char *use     = ap_main->get_aliname();
60    long        ali_len = GBT_get_alignment_len(GLOBAL_gb_main, use);
61    if (ali_len <= 1) {
62        aw_popup_exit("No valid alignment selected! Try again");
63    }
64
65    agt->tree_static->set_root_changed_callback(AWT_graphic_parsimony_root_changed, agt);
66}
67
68static double funktion_quadratisch(double wert, double *param_list, int param_anz) {
69    if (param_anz != 3) {
70        ap_assert(0); // wrong number of parameters
71        return 0;
72    }
73    return wert * wert * param_list[0] + wert * param_list[1] + param_list[2];
74}
75
76
77void ArbParsimony::kernighan_optimize_tree(AP_tree *at) {
78    GB_push_transaction(GLOBAL_gb_main);
79
80    long prevCombineCount = AP_sequence::combine_count();
81
82    AP_FLOAT       pars_start = get_root_node()->costs();
83    const AP_FLOAT pars_org   = pars_start;
84
85    int rek_deep_max = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/maxdepth");
86
87    AP_KL_FLAG funktype = (AP_KL_FLAG)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/function_type");
88
89    int param_anz;
90    double param_list[3];
91    double f_startx, f_maxy, f_maxx, f_max_deep;
92    f_max_deep = (double)rek_deep_max;                   //             x2
93    f_startx = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/start");
94    f_maxy = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/maxy");
95    f_maxx = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/maxx");
96
97    double (*funktion)(double wert, double *param_list, int param_anz);
98    switch (funktype) {
99        default:
100        case AP_QUADRAT_START:
101            funktion = funktion_quadratisch;
102            param_anz = 3;
103            param_list[2] = f_startx;
104            param_list[0] = (f_startx - f_maxy) / (f_maxx * f_maxx);
105            param_list[1] = -2.0 * param_list[0] * f_maxx;
106            break;
107        case AP_QUADRAT_MAX:    // parameter liste fuer quadratische gleichung (y =ax^2 +bx +c)
108            funktion = funktion_quadratisch;
109            param_anz = 3;
110            param_list[0] =  - f_maxy / ((f_max_deep -  f_maxx) * (f_max_deep - f_maxx));
111            param_list[1] =  -2.0 * param_list[0] * f_maxx;
112            param_list[2] =  f_maxy  + param_list[0] * f_maxx * f_maxx;
113            break;
114    }
115
116
117    AP_KL_FLAG searchflag=(AP_KL_FLAG)0;
118    if (*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/enable")) {
119        searchflag = AP_DYNAMIK;
120    }
121    if (*GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/enable")) {
122        searchflag = (AP_KL_FLAG)(searchflag|AP_STATIC);
123    }
124
125    int rek_breite[8];
126    rek_breite[0] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth0");
127    rek_breite[1] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth1");
128    rek_breite[2] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth2");
129    rek_breite[3] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth3");
130    rek_breite[4] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth4");
131    int rek_breite_anz = 5;
132
133    int       anzahl = (int)(*GBT_read_float(GLOBAL_gb_main, "genetic/kh/nodes")*at->count_leafs());
134    AP_tree **list   = at->getRandomNodes(anzahl);
135   
136    arb_progress progress(anzahl);
137
138    progress.subtitle(GBS_global_string("Old Parsimony: %.1f", pars_start));
139
140    GB_pop_transaction(GLOBAL_gb_main);
141
142    for (int i=0; i<anzahl && !progress.aborted(); i++) {
143        AP_tree_nlen *tree_elem = DOWNCAST(AP_tree_nlen*, list[i]); // @@@ pass 'at' as AP_tree_nlen
144
145        bool in_folded_group = tree_elem->gr.hidden ||
146            (tree_elem->father && tree_elem->get_father()->gr.hidden);
147
148        if (!in_folded_group) {
149            bool better_tree_found = false;
150            ap_main->push();
151            display_clear(funktion, param_list, param_anz, (int)pars_start, (int)rek_deep_max);
152
153            tree_elem->kernighan_rek(0,
154                                     rek_breite, rek_breite_anz, rek_deep_max,
155                                     funktion, param_list, param_anz,
156                                     pars_start,  pars_start, pars_org,
157                                     searchflag, &better_tree_found);
158
159            if (better_tree_found) {
160                ap_main->clear();
161                pars_start =  get_root_node()->costs();
162                progress.subtitle(GBS_global_string("New parsimony: %.1f (gain: %.1f)", pars_start, pars_org-pars_start));
163            }
164            else {
165                ap_main->pop();
166            }
167        }
168        progress.inc();
169    }
170    delete list;
171    printf("Combines: %li\n", AP_sequence::combine_count()-prevCombineCount);
172}
173
174
175
176void ArbParsimony::optimize_tree(AP_tree *at, arb_progress& progress) {
177    AP_tree        *oldrootleft  = get_root_node()->get_leftson();
178    AP_tree        *oldrootright = get_root_node()->get_rightson();
179    const AP_FLOAT  org_pars     = get_root_node()->costs();
180    AP_FLOAT        prev_pars    = org_pars;
181
182    progress.subtitle(GBS_global_string("Old parsimony: %.1f", org_pars));
183
184    while (!progress.aborted()) {
185        AP_FLOAT nni_pars = DOWNCAST(AP_tree_nlen*, at)->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false);
186
187        if (nni_pars == prev_pars) { // NNI did not reduce costs -> kern-lin
188            kernighan_optimize_tree(at);
189            AP_FLOAT ker_pars = get_root_node()->costs();
190            if (ker_pars == prev_pars) break; // kern-lin did not improve tree -> done
191            prev_pars = ker_pars;
192        }
193        else {
194            prev_pars = nni_pars;
195        }
196        progress.subtitle(GBS_global_string("New parsimony: %.1f (gain: %.1f)", prev_pars, org_pars-prev_pars));
197    }
198
199    if (oldrootleft->father == oldrootright) oldrootleft->set_root();
200    else oldrootright->set_root();
201
202    get_root_node()->costs();
203}
204
205AWT_graphic_parsimony::AWT_graphic_parsimony(ArbParsimony& parsimony_, GBDATA *gb_main_, AD_map_viewer_cb map_viewer_cb_)
206    : AWT_graphic_tree(parsimony_.get_awroot(), gb_main_, map_viewer_cb_),
207      parsimony(parsimony_)
208{}
209
210void ArbParsimony::generate_tree(WeightedFilter *pars_weighted_filter) {
211    AliView     *aliview   = pars_generate_aliview(pars_weighted_filter);
212    AP_sequence *seq_templ = 0;
213
214    GBDATA *gb_main = aliview->get_gb_main();
215    {
216        GB_transaction ta(gb_main);
217        bool           is_aa = GBT_is_alignment_protein(gb_main, aliview->get_aliname());
218
219        if (is_aa) seq_templ = new AP_sequence_protein(aliview);
220        else seq_templ       = new AP_sequence_parsimony(aliview);
221    }
222
223    tree = new AWT_graphic_parsimony(*this, aliview->get_gb_main(), PARS_map_viewer);
224    tree->init(new AP_TreeNlenNodeFactory, aliview, seq_templ, true, false);
225    ap_main->set_tree_root(tree);
226}
227
228AW_gc_manager AWT_graphic_parsimony::init_devices(AW_window *aww, AW_device *device, AWT_canvas* ntw) {
229    AW_init_color_group_defaults("arb_pars");
230
231    AW_gc_manager gc_manager =
232        AW_manage_GC(aww,
233                     ntw->get_gc_base_name(),
234                     device, AWT_GC_CURSOR, AWT_GC_MAX, /* AWT_GC_CURSOR+7, */ AW_GCM_DATA_AREA,
235                     makeWindowCallback(AWT_resize_cb, ntw),
236                     true,      // uses color groups
237                     "#AAAA55",
238
239                     // Important note :
240                     // Many gc indices are shared between ABR_NTREE and ARB_PARSIMONY
241                     // e.g. the tree drawing routines use same gc's for drawing both trees
242                     // (check AWT_dtree.cxx AWT_graphic_tree::init_devices)
243
244                     "Cursor$#FFFFFF",
245                     "Branch remarks$#DBE994",
246                     "+-Bootstrap$#DBE994",    "-B.(limited)$white",
247                     "--unused$#ff0000",
248                     "Marked$#FFFF00",
249                     "Some marked$#eeee88",
250                     "Not marked$black",
251                     "Zombies etc.$#cc5924",
252
253                     "--unused", "--unused", // these reserve the numbers which are used for probe colors in ARB_NTREE
254                     "--unused", "--unused", // (this is necessary because ARB_PARS and ARB_NTREE use the same tree painting routines)
255                     "--unused", "--unused",
256                     "--unused", "--unused",
257
258                     NULL);
259    return gc_manager;
260}
261
262void AWT_graphic_parsimony::show(AW_device *device) {
263    AP_tree_nlen *root_node = parsimony.get_root_node();
264    AW_awar      *awar_pars = aw_root->awar(AWAR_PARSIMONY);
265    AW_awar      *awar_best = aw_root->awar(AWAR_BEST_PARSIMONY);
266
267    long parsval = root_node ? root_node->costs() : 0;
268    awar_pars->write_int(parsval);
269
270    long best_pars = awar_best->read_int();
271    if (parsval < best_pars || 0==best_pars) awar_best->write_int(parsval);
272
273    AWT_graphic_tree::show(device);
274}
275
276void AWT_graphic_parsimony::handle_command(AW_device *device, AWT_graphic_event& event) {
277    ClickedTarget clicked(this, event.best_click());
278    bool          recalc_branchlengths_on_structure_change = true;
279
280    switch (event.cmd()) {
281        // @@@ something is designed completely wrong here!
282        // why do all commands close TA and reopen when done?
283
284        case AWT_MODE_NNI:
285            if (event.type()==AW_Mouse_Press) {
286                GB_pop_transaction(gb_main);
287                switch (event.button()) {
288                    case AW_BUTTON_LEFT: {
289                        if (clicked.node()) {
290                            arb_progress  progress("NNI optimize subtree");
291                            AP_tree_nlen *atn = DOWNCAST(AP_tree_nlen*, clicked.node());
292                            atn->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false);
293                            exports.save = 1;
294                            ASSERT_VALID_TREE(get_root_node());
295                        }
296                        break;
297                    }
298                    case AW_BUTTON_RIGHT: {
299                        arb_progress progress("NNI optimize tree");
300                        long         prevCombineCount = AP_sequence::combine_count();
301                       
302                        AP_tree_nlen *atn = DOWNCAST(AP_tree_nlen*, get_root_node());
303                        atn->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false);
304                        printf("Combines: %li\n", AP_sequence::combine_count()-prevCombineCount);
305
306                        exports.save = 1;
307                        ASSERT_VALID_TREE(get_root_node());
308                        break;
309                    }
310
311                    default: break;
312                }
313                GB_begin_transaction(gb_main);
314            }
315            break;
316        case AWT_MODE_KERNINGHAN:
317            if (event.type()==AW_Mouse_Press) {
318                GB_pop_transaction(gb_main);
319                switch (event.button()) {
320                    case AW_BUTTON_LEFT:
321                        if (clicked.node()) {
322                            arb_progress  progress("Kernighan-Lin optimize subtree");
323                            parsimony.kernighan_optimize_tree(clicked.node());
324                            this->exports.save = 1;
325                            ASSERT_VALID_TREE(get_root_node());
326                        }
327                        break;
328                    case AW_BUTTON_RIGHT: {
329                        arb_progress progress("Kernighan-Lin optimize tree");
330                        parsimony.kernighan_optimize_tree(get_root_node());
331                        this->exports.save = 1;
332                        ASSERT_VALID_TREE(get_root_node());
333                        break;
334                    }
335                    default: break;
336                }
337                GB_begin_transaction(gb_main);
338            }
339            break;
340        case AWT_MODE_OPTIMIZE:
341            if (event.type()==AW_Mouse_Press) {
342                GB_pop_transaction(gb_main);
343                switch (event.button()) {
344                    case AW_BUTTON_LEFT:
345                        if (clicked.node()) {
346                            arb_progress  progress("Optimizing subtree");
347                            parsimony.optimize_tree(clicked.node(), progress);
348                            this->exports.save = 1;
349                            ASSERT_VALID_TREE(get_root_node());
350                        }
351                        break;
352                    case AW_BUTTON_RIGHT: {
353                        arb_progress progress("Optimizing tree");
354                       
355                        parsimony.optimize_tree(get_root_node(), progress);
356                        this->exports.save = 1;
357                        ASSERT_VALID_TREE(get_root_node());
358                        break;
359                    }
360                    default: break;
361                }
362                GB_begin_transaction(gb_main);
363            }
364            break;
365
366        default:
367            recalc_branchlengths_on_structure_change = false;
368            // fall-through (modes listed below trigger branchlength calculation)
369        case AWT_MODE_MOVE:
370            AWT_graphic_tree::handle_command(device, event);
371            break;
372    }
373
374    if (exports.save == 1 && recalc_branchlengths_on_structure_change) {
375        arb_progress progress("Recalculating branch lengths");
376        rootEdge()->calc_branchlengths();
377        reorder_tree(BIG_BRANCHES_TO_TOP); // beautify after recalc_branch_lengths
378    }
379}
380
381
382// --------------------------------------------------------------------------------
383
384#ifdef UNIT_TESTS
385#include <arb_diff.h>
386#ifndef TEST_UNIT_H
387#include <test_unit.h>
388#endif
389
390void fake_AW_init_color_groups();
391
392struct fake_agt : public AWT_graphic_tree, virtual Noncopyable {
393    AP_sequence_parsimony *templ;
394
395    fake_agt()
396        : AWT_graphic_tree(NULL, GLOBAL_gb_main, NULL),
397          templ(NULL)
398    {
399    }
400    ~fake_agt() {
401        delete templ;
402    }
403    void init(AliView *aliview) {
404        fake_AW_init_color_groups(); // acts like no species has a color
405        delete templ;
406        templ = aliview->has_data() ? new AP_sequence_parsimony(aliview) : NULL;
407        AWT_graphic_tree::init(new AP_TreeNlenNodeFactory, aliview, templ, true, false);
408    }
409};
410
411class PARSIMONY_testenv : virtual Noncopyable {
412    GB_shell  shell;
413    AP_main   apMain;
414    fake_agt *agt;
415
416    void common_init(const char *dbname) {
417        GLOBAL_gb_main = NULL;
418        apMain.open(dbname);
419
420        TEST_EXPECT_NULL(ap_main);
421        ap_main = &apMain;
422
423        agt = new fake_agt;
424        apMain.set_tree_root(agt);
425    }
426
427public:
428    PARSIMONY_testenv(const char *dbname) {
429        common_init(dbname);
430        agt->init(new AliView(GLOBAL_gb_main));
431    }
432    PARSIMONY_testenv(const char *dbname, const char *aliName) {
433        common_init(dbname);
434        GB_transaction ta(GLOBAL_gb_main);
435        size_t aliLength = GBT_get_alignment_len(GLOBAL_gb_main, aliName);
436
437        AP_filter filter(aliLength);
438        if (!filter.is_invalid()) {
439            AP_weights weights(&filter);
440            agt->init(new AliView(GLOBAL_gb_main, filter, weights, aliName));
441        }
442    }
443    ~PARSIMONY_testenv() {
444        TEST_EXPECT_EQUAL(ap_main, &apMain);
445        ap_main = NULL;
446
447        delete agt;
448        GB_close(GLOBAL_gb_main);
449        GLOBAL_gb_main = NULL;
450    }
451
452    GB_ERROR load_tree(const char *tree_name) {
453        GB_transaction ta(GLOBAL_gb_main); // @@@ do inside AWT_graphic_tree::load?
454        return agt->load(GLOBAL_gb_main, tree_name, 0, 0);
455    }
456    AP_tree_nlen *tree_root() { return apMain.get_root_node(); }
457
458    void push() { apMain.push(); }
459    void pop() { apMain.pop(); }
460};
461
462void TEST_tree_modifications() {
463    PARSIMONY_testenv env("TEST_trees.arb");
464    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
465
466    {
467        AP_tree_nlen *root = env.tree_root();
468        TEST_REJECT_NULL(root);
469
470        AP_tree_edge::initialize(root);   // builds edges
471        TEST_ASSERT_VALID_TREE(root);
472
473        root->compute_tree();
474
475        // first check initial state:
476        {
477            AP_tree_members& root_info = root->gr;
478
479            TEST_EXPECT_EQUAL(root_info.grouped,             0);
480            TEST_EXPECT_EQUAL(root_info.hidden,              0);
481            TEST_EXPECT_EQUAL(root_info.has_marked_children, 1);
482            TEST_EXPECT_EQUAL(root_info.leaf_sum,            15);
483
484            TEST_EXPECT_SIMILAR(root_info.max_tree_depth, 1.624975, 0.000001);
485            TEST_EXPECT_SIMILAR(root_info.min_tree_depth, 0.341681, 0.000001);
486
487            GB_transaction ta(GLOBAL_gb_main);
488            GBT_mark_all(GLOBAL_gb_main, 0); // unmark all species
489            root->compute_tree();
490            TEST_EXPECT_EQUAL(root_info.has_marked_children, 0);
491        }
492
493
494#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"
495#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"
496#define B2_TOP "(((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010,CloPaste:0.179):0.131"
497#define B2_BOT "(CloPaste:0.179,(CloCarni:0.120,(CloButy2:0.009,CloButyr:0.000):0.564):0.010):0.131"
498
499
500#define B3_LEFT_TOP_SONS "(((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053,CelBiazo:0.059)"
501#define B3_TOP_SONS      B3_LEFT_TOP_SONS ":0.207,CytAquat:0.711"
502#define B3_TOP_SONS_CCR  "((CorAquat:0.187,CorGluta:0.522):0.053,CelBiazo:0.059):0.207,CytAquat:0.711" // CCR = CurCitre removed
503#define B3_TOP           "(" B3_TOP_SONS "):0.081"
504#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"
505
506
507        const char *top_topo    = "((" B1_TOP "," B2_TOP "):0.081," B3_TOP ");";
508        const char *edge_topo   = "((" B1_TOP "," B2_BOT "):0.081," B3_BOT ");";
509        const char *bottom_topo = "(" B3_BOT ",(" B2_BOT "," B1_BOT "):0.081);";
510
511        const char *radial_topo  =
512            "(((CloPaste:0.179,((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010):0.131,"
513            "((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,"
514            "((CelBiazo:0.059,((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053):0.207,CytAquat:0.711):0.081);";
515        const char *radial_topo2 =
516            "(((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,"
517            "(CytAquat:0.711," B3_LEFT_TOP_SONS ":0.207):0.081);";
518
519        // expect that no mode reproduces another mode:
520        TEST_EXPECT_DIFFERENT(top_topo,    edge_topo);
521        TEST_EXPECT_DIFFERENT(top_topo,    bottom_topo);
522        TEST_EXPECT_DIFFERENT(top_topo,    radial_topo);
523        TEST_EXPECT_DIFFERENT(top_topo,    radial_topo2);
524        TEST_EXPECT_DIFFERENT(edge_topo,   bottom_topo);
525        TEST_EXPECT_DIFFERENT(edge_topo,   radial_topo);
526        TEST_EXPECT_DIFFERENT(edge_topo,   radial_topo2);
527        TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo);
528        TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo2);
529        TEST_EXPECT_DIFFERENT(radial_topo, radial_topo2);
530
531        env.push(); // 1st stack level (=top_topo)
532
533        TEST_ASSERT_VALID_TREE(root);
534
535        TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
536        // test reorder_tree:
537        root->reorder_tree(BIG_BRANCHES_TO_EDGE);     TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo);    env.push(); // 2nd stack level (=edge_topo)
538        root->reorder_tree(BIG_BRANCHES_TO_BOTTOM);   TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo);  env.push(); // 3rd stack level (=bottom_topo)
539        root->reorder_tree(BIG_BRANCHES_TO_CENTER);   TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo);
540        root->reorder_tree(BIG_BRANCHES_ALTERNATING); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo2);
541        root->reorder_tree(BIG_BRANCHES_TO_TOP);      TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
542
543        TEST_ASSERT_VALID_TREE(root);
544
545        // test set root:
546        AP_tree_nlen *CloTyrob = root->findLeafNamed("CloTyrob");
547        TEST_REJECT_NULL(CloTyrob);
548
549        ARB_edge rootEdge(root->get_leftson(), root->get_rightson());
550        CloTyrob->set_root();
551
552        TEST_ASSERT_VALID_TREE(root);
553
554        const char *rootAtCloTyrob_topo =
555            "(CloTyrob:0.004,"
556            "(((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,"
557            "((((" B3_TOP_SONS "):0.162," B2_TOP "):0.124,CloBifer:0.388):0.057,CloInnoc:0.371):0.274):0.004);";
558
559        TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo);
560        env.push(); // 4th stack level (=rootAtCloTyrob_topo)
561
562        TEST_ASSERT_VALID_TREE(root);
563
564        AP_tree_nlen *CelBiazoFather = root->findLeafNamed("CelBiazo")->get_father();
565        TEST_REJECT_NULL(CelBiazoFather);
566        CelBiazoFather->set_root();
567
568        const char *rootAtCelBiazoFather_topo = "(" B3_LEFT_TOP_SONS ":0.104,((" B1_TOP "," B2_TOP "):0.162,CytAquat:0.711):0.104);";
569        TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCelBiazoFather_topo);
570
571        TEST_ASSERT_VALID_TREE(root);
572
573        ARB_edge oldRootEdge(rootEdge.source(), rootEdge.dest());
574        DOWNCAST(AP_tree_nlen*,oldRootEdge.son())->set_root();
575
576        const char *rootSetBack_topo = top_topo;
577        TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo);
578        env.push(); // 5th stack level (=rootSetBack_topo)
579
580        TEST_ASSERT_VALID_TREE(root);
581
582        // test remove:
583        AP_tree_nlen *CurCitre = root->findLeafNamed("CurCitre");
584        TEST_REJECT_NULL(CurCitre);
585        TEST_REJECT_NULL(CurCitre->get_father());
586
587        CurCitre->remove();
588        const char *CurCitre_removed_topo = "((" B1_TOP "," B2_TOP "):0.081,(" B3_TOP_SONS_CCR "):0.081);";
589        // ------------------------------------------------------------------- ^^^ = B3_TOP_SONS minus CurCitre
590        TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo);
591
592        TEST_ASSERT_VALID_TREE(root);
593        TEST_ASSERT_VALID_TREE(CurCitre);
594
595        TEST_EXPECT_EQUAL(root->gr.leaf_sum, 15); // out of date
596        root->compute_tree();
597        TEST_EXPECT_EQUAL(root->gr.leaf_sum, 14);
598
599        env.push(); // 6th stack level (=CurCitre_removed_topo)
600
601        TEST_ASSERT_VALID_TREE(root);
602
603        // test insert:
604        AP_tree_nlen *CloCarni = root->findLeafNamed("CloCarni");
605        TEST_REJECT_NULL(CloCarni);
606        CurCitre->insert(CloCarni); // this creates two extra edges (not destroyed by destroy() below) and one extra node
607
608        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);";
609        TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_inserted_topo);
610
611        AP_tree_nlen *node_del_manually  = CurCitre->get_father();
612        AP_tree_edge *edge1_del_manually = CurCitre->edgeTo(node_del_manually);
613        AP_tree_edge *edge2_del_manually = CurCitre->get_brother()->edgeTo(node_del_manually);
614
615        TEST_ASSERT_VALID_TREE(root);
616
617        // now check pops:
618        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo);
619        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo);
620        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo);
621        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo);
622        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo);
623        env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
624
625        TEST_ASSERT_VALID_TREE(root);
626
627        AP_tree_edge::destroy(root);
628
629        // delete memory allocated by insert() above and lost due to pop()s
630        delete edge1_del_manually;
631        delete edge2_del_manually;
632
633        node_del_manually->forget_origin();
634        node_del_manually->father   = NULL;
635        node_del_manually->leftson  = NULL;
636        node_del_manually->rightson = NULL;
637        delete node_del_manually;
638    }
639}
640
641// @@@ Tests wanted:
642// - NNI
643// - tree optimize
644// - ...
645
646void TEST_calc_bootstraps() {
647    PARSIMONY_testenv env("TEST_trees.arb", "ali_5s");
648    TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
649
650    const char *bs_origi_topo = "(((((((CloTyro3,CloTyro4)'40%',CloTyro2)'0%',CloTyrob)'97%',CloInnoc)'0%',CloBifer)'53%',(((CloButy2,CloButyr)'100%',CloCarni)'33%',CloPaste)'97%')'100%',((((CorAquat,CurCitre)'100%',CorGluta)'17%',CelBiazo)'40%',CytAquat)'100%');";
651    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%');";
652    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%');";
653
654    {
655        AP_tree_nlen *root = env.tree_root();
656        TEST_REJECT_NULL(root);
657
658        AP_tree_edge::initialize(root);   // builds edges
659        TEST_EXPECT_EQUAL(root, rootNode()); // need tree-access via global 'ap_main' (too much code is based on that)
660
661        AP_tree_edge *root_edge = rootEdge();
662        TEST_REJECT_NULL(root_edge);
663
664        root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_origi_topo);
665
666        root_edge->nni_rek(-1, false, AP_BL_MODE(AP_BL_BL_ONLY|AP_BL_BOOTSTRAP_LIMIT),    NULL); root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_limit_topo);
667        root_edge->nni_rek(-1, false, AP_BL_MODE(AP_BL_BL_ONLY|AP_BL_BOOTSTRAP_ESTIMATE), NULL); root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_estim_topo);
668
669        TEST_EXPECT_EQUAL(env.tree_root(), root);
670        AP_tree_edge::destroy(root);
671    }
672}
673
674void TEST_tree_remove_add_all() {
675    // reproduces crash as described in #527
676    PARSIMONY_testenv env("TEST_trees.arb", "ali_5s");
677    TEST_EXPECT_NO_ERROR(env.load_tree("tree_nj"));
678
679    const int     LEAFS     = 6;
680    AP_tree_nlen *leaf[LEAFS];
681    const char *name[LEAFS] = {
682        "CloButy2",
683        "CloButyr",
684        "CytAquat",
685        "CorAquat",
686        "CurCitre",
687        "CorGluta",
688    };
689
690    AP_tree_nlen *root = env.tree_root();
691    TEST_REJECT_NULL(root);
692
693    for (int i = 0; i<LEAFS; ++i) {
694        leaf[i] = root->findLeafNamed(name[i]);
695        TEST_REJECT_NULL(leaf[i]);
696    }
697
698    AP_tree_edge::initialize(root);   // builds edges
699    TEST_EXPECT_EQUAL(root, rootNode()); // need tree-access via global 'ap_main' (too much code is based on that)
700
701    AP_tree_root *troot = leaf[0]->get_tree_root();
702    TEST_REJECT_NULL(troot);
703
704    // Note: following loop leaks father nodes and edges
705    for (int i = 0; i<LEAFS-1; ++i) { // removing the second to last leaf, "removes" both remaining leafs
706        TEST_ASSERT_VALID_TREE(root);
707        leaf[i]->remove();
708        TEST_ASSERT_VALID_TREE(leaf[i]);
709    }
710    leaf[LEAFS-1]->father = NULL; // correct final leaf (not removed regularily)
711
712    leaf[0]->initial_insert(leaf[1], troot);
713    for (int i = 2; i<LEAFS; ++i) {
714        TEST_ASSERT_VALID_TREE(leaf[i-1]);
715        TEST_ASSERT_VALID_TREE(leaf[i]);
716        leaf[i]->insert(leaf[i-1]);
717    }
718}
719
720#endif // UNIT_TESTS
721
722// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.