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