source: tags/ms_r17q3/ARBDB/TreeNode.cxx

Last change on this file was 15823, checked in by westram, 8 years ago
  • partial merge of [15807] from 'group' into 'trunk'
  • 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      : TreeNode.cxx                                       //
4//   Purpose   :                                                    //
5//                                                                  //
6//   Coded by Ralf Westram (coder@reallysoft.de) in December 2013   //
7//   Institute of Microbiology (Technical University Munich)        //
8//   http://www.arb-home.de/                                        //
9//                                                                  //
10// ================================================================ //
11
12#include "TreeNode.h"
13#include <arb_progress.h>
14#include <map>
15#include <set>
16#include <cmath> // needed with 4.4.3 (but not with 4.7.3)
17
18// ------------------
19//      TreeRoot
20
21TreeRoot::~TreeRoot() {
22    deleteWithNodes = false; // avoid recursive call of ~TreeRoot (obsolete?)
23    rt_assert(!rootNode);    // you have to call TreeRoot::predelete() before dtor! you can do this is dtor of that derived class, which defines makeNode/destroyNode
24    // Note: destroying nodes from here is impossible (leads to pure virtual call, as derived class instance of 'this' is already under destruction)
25}
26
27void TreeRoot::change_root(TreeNode *oldroot, TreeNode *newroot) {
28    rt_assert(rootNode == oldroot);
29    rt_assert(implicated(newroot, !newroot->father));
30    rootNode = newroot;
31
32    if (oldroot && oldroot->get_tree_root() && !oldroot->is_inside(newroot)) oldroot->set_tree_root(0); // unlink from this
33    if (newroot && newroot->get_tree_root() != this) newroot->set_tree_root(this); // link to this
34}
35
36// --------------------
37//      TreeNode
38
39#if defined(PROVIDE_TREE_STRUCTURE_TESTS)
40
41Validity TreeNode::is_valid() const {
42    rt_assert(knownNonNull(this));
43    Validity valid;
44
45    TreeRoot *troot = get_tree_root();
46    if (troot) {
47        if (is_leaf) {
48            valid = Validity(!rightson && !leftson, "leaf has son");
49        }
50        else {
51            valid            = Validity(rightson && leftson, "inner node lacks sons");
52            if (valid) valid = get_rightson()->is_valid();
53            if (valid) valid = get_leftson()->is_valid();
54        }
55        if (father) {
56            if (valid) valid = Validity(is_inside(get_father()), "node not inside father subtree");
57            if (valid) valid = Validity(troot->get_root_node()->is_anchestor_of(this), "root is not nodes anchestor");
58            if (valid) valid = Validity(get_father()->get_tree_root() == troot, "node and father have different TreeRoot");
59        }
60        else {
61            if (valid) valid = Validity(troot->get_root_node() == this, "node has no father, but isn't root-node");
62            if (valid) valid = Validity(!is_leaf, "root-node is leaf"); // leaf@root (tree has to have at least 2 leafs)
63        }
64    }
65    else { // removed node (may be incomplete)
66        if (is_leaf) {
67            valid = Validity(!rightson && !leftson, "(removed) leaf has son");
68        }
69        else {
70            if (rightson)         valid = get_rightson()->is_valid();
71            if (leftson && valid) valid = get_leftson()->is_valid();
72        }
73        if (father) {
74            if (valid) valid = Validity(is_inside(get_father()), "(removed) node not inside father subtree");
75            if (valid) valid = Validity(get_father()->get_tree_root() == troot, "(removed) node and father have different TreeRoot");
76        }
77    }
78    return valid;
79}
80#endif // PROVIDE_TREE_STRUCTURE_TESTS
81
82void TreeNode::set_tree_root(TreeRoot *new_root) {
83    if (tree_root != new_root) {
84        tree_root = new_root;
85        if (leftson) get_leftson()->set_tree_root(new_root);
86        if (rightson) get_rightson()->set_tree_root(new_root);
87    }
88}
89
90void TreeNode::reorder_subtree(TreeOrder mode) {
91    static const char *smallest_leafname; // has to be set to the alphabetically smallest name (when function exits)
92
93    if (is_leaf) {
94        smallest_leafname = name;
95    }
96    else {
97        int leftsize  = get_leftson() ->get_leaf_count();
98        int rightsize = get_rightson()->get_leaf_count();
99
100        {
101            bool big_at_top    = leftsize>rightsize;
102            bool big_at_bottom = leftsize<rightsize;
103            bool swap_branches = (mode&ORDER_BIG_DOWN) ? big_at_top : big_at_bottom;
104            if (swap_branches) swap_sons();
105        }
106
107        TreeOrder lmode, rmode;
108        if (mode & (ORDER_BIG_TO_EDGE|ORDER_BIG_TO_CENTER)) { // symmetric
109            if (mode & ORDER_ALTERNATING) mode = TreeOrder(mode^(ORDER_BIG_TO_EDGE|ORDER_BIG_TO_CENTER));
110
111            if (mode & ORDER_BIG_TO_CENTER) {
112                lmode = TreeOrder(mode | ORDER_BIG_DOWN);
113                rmode = TreeOrder(mode & ~ORDER_BIG_DOWN);
114            }
115            else {
116                lmode = TreeOrder(mode & ~ORDER_BIG_DOWN);
117                rmode = TreeOrder(mode | ORDER_BIG_DOWN);
118            }
119        }
120        else { // asymmetric
121            if (mode & ORDER_ALTERNATING) mode = TreeOrder(mode^ORDER_BIG_DOWN);
122
123            lmode = mode;
124            rmode = mode;
125        }
126
127        get_leftson()->reorder_subtree(lmode);
128        const char *leftleafname = smallest_leafname;
129
130        get_rightson()->reorder_subtree(rmode);
131        const char *rightleafname = smallest_leafname;
132
133        if (leftleafname && rightleafname) {
134            int name_cmp = strcmp(leftleafname, rightleafname);
135            if (name_cmp <= 0) {
136                smallest_leafname = leftleafname;
137            }
138            else {
139                smallest_leafname = rightleafname;
140                if (leftsize == rightsize) { // if sizes of subtrees are equal and rightleafname<leftleafname -> swap branches
141                    const char *smallest_leafname_save = smallest_leafname;
142
143                    swap_sons();
144                    get_leftson ()->reorder_subtree(lmode); rt_assert(strcmp(smallest_leafname, rightleafname)== 0);
145                    get_rightson()->reorder_subtree(rmode); rt_assert(strcmp(smallest_leafname, leftleafname) == 0);
146
147                    smallest_leafname = smallest_leafname_save;
148                }
149            }
150        }
151    }
152    rt_assert(smallest_leafname);
153}
154
155void TreeNode::reorder_tree(TreeOrder mode) {
156    /*! beautify tree (does not change topology, only swaps branches)
157     */
158    compute_tree();
159    reorder_subtree(mode);
160}
161
162void TreeNode::rotate_subtree() {
163    if (!is_leaf) {
164        swap_sons();
165        get_leftson()->rotate_subtree();
166        get_rightson()->rotate_subtree();
167    }
168}
169
170void TreeNode::set_root() {
171    /*! set the root at parent edge of this
172     * pointers to tree-nodes remain valid, but all parent-nodes of this change their meaning
173     * (afterwards they will point to [father+brother] instead of [this+brother])
174     * esp. pointers to the root-node will still point to the root-node (which only changed children)
175     */
176
177    if (at_root()) return; // already root
178
179    TreeNode *old_root    = get_root_node();
180    TreeNode *old_brother = is_inside(old_root->get_leftson()) ? old_root->get_rightson() : old_root->get_leftson();
181
182    // move remark branches to top
183    {
184        char *remark = nulldup(get_remark());
185        for (TreeNode *node = this; node->father; node = node->get_father()) {
186            remark = node->swap_remark(remark);
187        }
188        free(remark);
189    }
190
191    GBT_LEN old_root_len = old_root->leftlen + old_root->rightlen;
192
193    // new node & this init
194    old_root->leftson  = this;
195    old_root->rightson = father; // will be set later
196
197    if (father->leftson == this) {
198        old_root->leftlen = old_root->rightlen = father->leftlen*.5;
199    }
200    else {
201        old_root->leftlen = old_root->rightlen = father->rightlen*.5;
202    }
203
204    TreeNode *next = get_father()->get_father();
205    TreeNode *prev = old_root;
206    TreeNode *pntr = get_father();
207
208    if (father->leftson == this)    father->leftson = old_root; // to set the flag correctly
209
210    // loop from father to son of root, rotate tree
211    while  (next->father) {
212        double len = (next->leftson == pntr) ? next->leftlen : next->rightlen;
213
214        if (pntr->leftson == prev) {
215            pntr->leftson = next;
216            pntr->leftlen = len;
217        }
218        else {
219            pntr->rightson = next;
220            pntr->rightlen = len;
221        }
222
223        pntr->father = prev;
224        prev         = pntr;
225        pntr         = next;
226        next         = next->get_father();
227    }
228    // now next points to the old root, which has been destroyed above
229    //
230    // pointer at oldroot
231    // pntr == brother before old root == next
232
233    if (pntr->leftson == prev) {
234        pntr->leftlen = old_root_len;
235        pntr->leftson = old_brother;
236    }
237    else {
238        pntr->rightlen = old_root_len;
239        pntr->rightson = old_brother;
240    }
241
242    old_brother->father = pntr;
243    pntr->father        = prev;
244    father              = old_root;
245
246    rt_assert(get_root_node() == old_root);
247}
248
249TreeNode *TreeNode::findLeafNamed(const char *wantedName) {
250    TreeNode *found = NULL;
251    if (is_leaf) {
252        if (name && strcmp(name, wantedName) == 0) found = this;
253    }
254    else {
255        found             = get_leftson()->findLeafNamed(wantedName);
256        if (!found) found = get_rightson()->findLeafNamed(wantedName);
257    }
258    return found;
259}
260
261// ----------------------------
262//      find_innermost_edge
263
264class NodeLeafDistance {
265    GBT_LEN downdist, updist;
266    enum { NLD_NODIST = 0, NLD_DOWNDIST, NLD_BOTHDIST } state;
267
268public:
269
270    NodeLeafDistance()
271        : downdist(-1.0),
272          updist(-1.0),
273          state(NLD_NODIST)
274    {}
275
276    GBT_LEN get_downdist() const { rt_assert(state >= NLD_DOWNDIST); return downdist; }
277    void set_downdist(GBT_LEN DownDist) {
278        if (state < NLD_DOWNDIST) state = NLD_DOWNDIST;
279        downdist = DownDist;
280    }
281
282    GBT_LEN get_updist() const { rt_assert(state >= NLD_BOTHDIST); return updist; }
283    void set_updist(GBT_LEN UpDist) {
284        if (state < NLD_BOTHDIST) state = NLD_BOTHDIST;
285        updist = UpDist;
286    }
287
288};
289
290class EdgeFinder {
291    std::map<TreeNode*, NodeLeafDistance> data; // maximum distance to farthest leaf
292
293    ARB_edge innermost;
294    double   min_distdiff; // abs diff between up- and downdiff
295
296    GBT_LEN calc_distdiff(GBT_LEN d1, GBT_LEN d2) { return fabs(d1-d2); }
297
298    void insert_tree(TreeNode *node) {
299        if (node->is_leaf) {
300            data[node].set_downdist(0.0);
301        }
302        else {
303            insert_tree(node->get_leftson());
304            insert_tree(node->get_rightson());
305
306            data[node].set_downdist(std::max(data[node->get_leftson()].get_downdist()+node->leftlen,
307                                             data[node->get_rightson()].get_downdist()+node->rightlen));
308        }
309    }
310
311    void findBetterEdge_sub(TreeNode *node) {
312        TreeNode *father  = node->get_father();
313        TreeNode *brother = node->get_brother();
314
315        GBT_LEN len      = node->get_branchlength();
316        GBT_LEN brothLen = brother->get_branchlength();
317
318        GBT_LEN upDist   = std::max(data[father].get_updist(), data[brother].get_downdist()+brothLen);
319        GBT_LEN downDist = data[node].get_downdist();
320
321        {
322            GBT_LEN edge_dd = calc_distdiff(upDist, downDist);
323            if (edge_dd<min_distdiff) { // found better edge
324                innermost    = ARB_edge(node, father);
325                min_distdiff = edge_dd;
326            }
327        }
328
329        data[node].set_updist(upDist+len);
330
331        if (!node->is_leaf) {
332            findBetterEdge_sub(node->get_leftson());
333            findBetterEdge_sub(node->get_rightson());
334        }
335    }
336
337    void findBetterEdge(TreeNode *node) {
338        if (!node->is_leaf) {
339            findBetterEdge_sub(node->get_leftson());
340            findBetterEdge_sub(node->get_rightson());
341        }
342    }
343
344public:
345    EdgeFinder(TreeNode *rootNode)
346        : innermost(rootNode->get_leftson(), rootNode->get_rightson()) // root-edge
347    {
348        insert_tree(rootNode);
349
350        TreeNode *lson = rootNode->get_leftson();
351        TreeNode *rson = rootNode->get_rightson();
352
353        GBT_LEN rootEdgeLen = rootNode->leftlen + rootNode->rightlen;
354
355        GBT_LEN lddist = data[lson].get_downdist();
356        GBT_LEN rddist = data[rson].get_downdist();
357
358        data[lson].set_updist(rddist+rootEdgeLen);
359        data[rson].set_updist(lddist+rootEdgeLen);
360
361        min_distdiff = calc_distdiff(lddist, rddist);
362
363        findBetterEdge(lson);
364        findBetterEdge(rson);
365    }
366
367    const ARB_edge& innermost_edge() const { return innermost; }
368};
369
370ARB_edge TreeRoot::find_innermost_edge() {
371    EdgeFinder edgeFinder(get_root_node());
372    return edgeFinder.innermost_edge();
373}
374
375// ------------------------
376//      multifurcation
377
378class TreeNode::LengthCollector {
379    typedef std::map<TreeNode*,GBT_LEN> LengthMap;
380    typedef std::set<TreeNode*>         NodeSet;
381
382    LengthMap eliminatedParentLength;
383    LengthMap addedParentLength;
384
385public:
386    void eliminate_parent_edge(TreeNode *node) {
387        rt_assert(!node->is_root_node());
388        eliminatedParentLength[node] += parentEdge(node).eliminate();
389    }
390
391    void add_parent_length(TreeNode *node, GBT_LEN addLen) {
392        rt_assert(!node->is_root_node());
393        addedParentLength[node] += addLen;
394    }
395
396    void independent_distribution(bool useProgress) {
397        // step 2: (see caller)
398        arb_progress *progress    = NULL;
399        int           redistCount = 0;
400        if (useProgress) progress = new arb_progress("Distributing eliminated lengths", eliminatedParentLength.size());
401
402        while (!eliminatedParentLength.empty()) { // got eliminated lengths which need to be distributed
403            for (LengthMap::iterator from = eliminatedParentLength.begin(); from != eliminatedParentLength.end(); ++from) {
404                ARB_edge elimEdge = parentEdge(from->first);
405                GBT_LEN  elimLen  = from->second;
406
407                elimEdge.virtually_distribute_length(elimLen, *this);
408                if (progress) ++*progress;
409            }
410            eliminatedParentLength.clear(); // all distributed!
411
412            // handle special cases where distributed length is negative and results in negative destination branches.
413            // Avoid generating negative dest. branchlengths by
414            // - eliminating the dest. branch
415            // - redistributing the additional (negative) length (may cause additional negative lengths on other dest. branches)
416
417            NodeSet handled;
418            for (LengthMap::iterator to = addedParentLength.begin(); to != addedParentLength.end(); ++to) {
419                ARB_edge affectedEdge     = parentEdge(to->first);
420                GBT_LEN  additionalLen    = to->second;
421                double   effective_length = affectedEdge.length() + additionalLen;
422
423                if (effective_length<=0.0) { // negative or zero
424                    affectedEdge.set_length(effective_length);
425                    eliminate_parent_edge(to->first); // adds entry to eliminatedParentLength and causes another additional loop
426                    handled.insert(to->first);
427                }
428            }
429
430            if (progress && !eliminatedParentLength.empty()) {
431                delete progress;
432                ++redistCount;
433                progress = new arb_progress(GBS_global_string("Redistributing negative lengths (#%i)", redistCount), eliminatedParentLength.size());
434            }
435
436            // remove all redistributed nodes
437            for (NodeSet::iterator del = handled.begin(); del != handled.end(); ++del) {
438                addedParentLength.erase(*del);
439            }
440        }
441
442        // step 3:
443        for (LengthMap::iterator to = addedParentLength.begin(); to != addedParentLength.end(); ++to) {
444            ARB_edge affectedEdge     = parentEdge(to->first);
445            GBT_LEN  additionalLen    = to->second;
446            double   effective_length = affectedEdge.length() + additionalLen;
447
448            affectedEdge.set_length(effective_length);
449        }
450
451        if (progress) delete progress;
452    }
453};
454
455GBT_LEN ARB_edge::adjacent_distance() const {
456    //! return length of edges starting from this->dest()
457
458    if (is_edge_to_leaf()) return 0.0;
459    return next().length_or_adjacent_distance() + counter_next().length_or_adjacent_distance();
460}
461
462void ARB_edge::virtually_add_or_distribute_length_forward(GBT_LEN len, TreeNode::LengthCollector& collect) const {
463    rt_assert(!is_nan_or_inf(len));
464    if (length() > 0.0) {
465        collect.add_parent_length(son(), len);
466    }
467    else {
468        if (len != 0.0) virtually_distribute_length_forward(len, collect);
469    }
470}
471
472
473void ARB_edge::virtually_distribute_length_forward(GBT_LEN len, TreeNode::LengthCollector& collect) const {
474    /*! distribute length to edges adjacent in edge direction (i.e. edges starting from this->dest())
475     * Split 'len' proportional to adjacent edges lengths.
476     *
477     * Note: length will not be distributed to tree-struction itself (yet), but collected in 'collect'
478     */
479
480    rt_assert(is_normal(len));
481    rt_assert(!is_edge_to_leaf()); // cannot forward anything (nothing beyond leafs)
482
483    ARB_edge e1 = next();
484    ARB_edge e2 = counter_next();
485
486    GBT_LEN d1 = e1.length_or_adjacent_distance();
487    GBT_LEN d2 = e2.length_or_adjacent_distance();
488
489    len = len/(d1+d2);
490
491    e1.virtually_add_or_distribute_length_forward(len*d1, collect);
492    e2.virtually_add_or_distribute_length_forward(len*d2, collect);
493}
494
495void ARB_edge::virtually_distribute_length(GBT_LEN len, TreeNode::LengthCollector& collect) const {
496    /*! distribute length to all adjacent edges.
497     * Longer edges receive more than shorter ones.
498     *
499     * Edges with length zero will not be changed, instead both edges "beyond"
500     * the edge will be affected (they will be affected equally to direct edges,
501     * because edges at multifurcations are considered to BE direct edges).
502     *
503     * Note: length will not be distributed to tree-struction itself (yet), but collected in 'collect'
504     */
505
506    ARB_edge backEdge = inverse();
507    GBT_LEN len_fwd, len_bwd;
508    {
509        GBT_LEN dist_fwd = adjacent_distance();
510        GBT_LEN dist_bwd = backEdge.adjacent_distance();
511        GBT_LEN lenW     = len/(dist_fwd+dist_bwd);
512        len_fwd          = lenW*dist_fwd;
513        len_bwd          = lenW*dist_bwd;
514
515    }
516
517    if (is_normal(len_fwd)) virtually_distribute_length_forward(len_fwd, collect);
518    if (is_normal(len_bwd)) backEdge.virtually_distribute_length_forward(len_bwd, collect);
519}
520
521void TreeNode::eliminate_and_collect(const multifurc_limits& below, LengthCollector& collect) {
522    /*! eliminate edges specified by 'below' and
523     * store their length in 'collect' for later distribution.
524     */
525    rt_assert(!is_root_node());
526    if (!is_leaf || below.applyAtLeafs) {
527        double value;
528        switch (parse_bootstrap(value)) {
529            case REMARK_NONE:
530                value = 100.0;
531                // fall-through
532            case REMARK_BOOTSTRAP:
533                if (value<below.bootstrap && get_branchlength_unrooted()<below.branchlength) {
534                    collect.eliminate_parent_edge(this);
535                }
536                break;
537
538            case REMARK_OTHER: break;
539        }
540    }
541
542    if (!is_leaf) {
543        get_leftson() ->eliminate_and_collect(below, collect);
544        get_rightson()->eliminate_and_collect(below, collect);
545    }
546}
547
548void ARB_edge::multifurcate() {
549    /*! eliminate edge and distribute length to adjacent edges
550     * - sets its length to zero
551     * - removes remark (bootstrap)
552     * - distributes length to neighbour-branches
553     */
554    TreeNode::LengthCollector collector;
555    collector.eliminate_parent_edge(son());
556    collector.independent_distribution(false);
557}
558void TreeNode::multifurcate() {
559    /*! eliminate branch from 'this' to 'father' (or brother @ root)
560     * @see ARB_edge::multifurcate()
561     */
562    rt_assert(father); // cannot multifurcate at root; call with son of root to multifurcate root-edge
563    if (father) parentEdge(this).multifurcate();
564}
565
566void TreeNode::set_branchlength_preserving(GBT_LEN new_len) {
567    /*! set branchlength to 'new_len' while preserving overall distance in tree.
568     *
569     * Always works on unrooted tree (i.e. lengths @ root are treated correctly).
570     * Length is preserved as in multifurcate()
571     */
572
573    GBT_LEN old_len = get_branchlength_unrooted();
574    GBT_LEN change  = new_len-old_len;
575
576    char *old_remark = nulldup(get_remark());
577
578    // distribute the negative 'change' to neighbours:
579    set_branchlength_unrooted(-change);
580    multifurcate();
581
582    set_branchlength_unrooted(new_len);
583    use_as_remark(old_remark); // restore remark (was removed by multifurcate())
584}
585
586void TreeNode::multifurcate_whole_tree(const multifurc_limits& below) {
587    /*! multifurcate all branches specified by 'below'
588     * - step 1: eliminate all branches, store eliminated lengths
589     * - step 2: calculate length distribution for all adjacent branches (proportionally to original length of each branch)
590     * - step 3: apply distributed length
591     */
592    TreeNode        *root = get_root_node();
593    LengthCollector  collector;
594    arb_progress progress("Multifurcating tree");
595
596    // step 1:
597    progress.subtitle("Collecting branches to eliminate");
598    root->get_leftson()->eliminate_and_collect(below, collector);
599    root->get_rightson()->eliminate_and_collect(below, collector);
600    // root-edge is handled twice by the above calls.
601    // Unproblematic: first call will eliminate root-edge, so second call will do nothing.
602
603    // step 2 and 3:
604    collector.independent_distribution(true);
605}
606
607TreeNode::bs100_mode TreeNode::toggle_bootstrap100(bs100_mode mode) {
608    if (!is_leaf) {
609        if (!is_root_node()) {
610            double bootstrap;
611            switch (parse_bootstrap(bootstrap)) {
612                case REMARK_NONE:
613                case REMARK_OTHER:
614                    switch (mode) {
615                        case BS_UNDECIDED: mode = BS_INSERT;
616                        case BS_INSERT: set_bootstrap(100);
617                        case BS_REMOVE: break;
618                    }
619                    break;
620                case REMARK_BOOTSTRAP:
621                    if (bootstrap >= 99.5) {
622                        switch (mode) {
623                            case BS_UNDECIDED: mode = BS_REMOVE;
624                            case BS_REMOVE: remove_remark();
625                            case BS_INSERT: break;
626                        }
627                    }
628                    break;
629            }
630        }
631
632        mode = get_leftson()->toggle_bootstrap100(mode);
633        mode = get_rightson()->toggle_bootstrap100(mode);
634    }
635    return mode;
636}
637void TreeNode::remove_bootstrap() {
638    freenull(remark_branch);
639    if (!is_leaf) {
640        get_leftson()->remove_bootstrap();
641        get_rightson()->remove_bootstrap();
642    }
643}
644void TreeNode::reset_branchlengths() {
645    if (!is_leaf) {
646        leftlen = rightlen = DEFAULT_BRANCH_LENGTH;
647
648        get_leftson()->reset_branchlengths();
649        get_rightson()->reset_branchlengths();
650    }
651}
652
653void TreeNode::scale_branchlengths(double factor) {
654    if (!is_leaf) {
655        leftlen  *= factor;
656        rightlen *= factor;
657
658        get_leftson()->scale_branchlengths(factor);
659        get_rightson()->scale_branchlengths(factor);
660    }
661}
662
663GBT_LEN TreeNode::sum_child_lengths() const {
664    if (is_leaf) return 0.0;
665    return
666        leftlen +
667        rightlen +
668        get_leftson()->sum_child_lengths() +
669        get_rightson()->sum_child_lengths();
670}
671
672void TreeNode::bootstrap2branchlen() {
673    //! copy bootstraps to branchlengths
674    if (is_leaf) {
675        set_branchlength_unrooted(DEFAULT_BRANCH_LENGTH);
676    }
677    else {
678        if (father) {
679            double         bootstrap;
680            GBT_RemarkType rtype = parse_bootstrap(bootstrap);
681
682            if (rtype == REMARK_BOOTSTRAP) {
683                double len = bootstrap/100.0;
684                set_branchlength_unrooted(len);
685            }
686            else {
687                set_branchlength_unrooted(1.0); // no bootstrap means "100%"
688            }
689        }
690        get_leftson()->bootstrap2branchlen();
691        get_rightson()->bootstrap2branchlen();
692    }
693}
694
695void TreeNode::branchlen2bootstrap() {
696    //! copy branchlengths to bootstraps
697    remove_remark();
698    if (!is_leaf) {
699        if (!is_root_node()) {
700            set_bootstrap(get_branchlength_unrooted()*100.0);
701        }
702        get_leftson()->branchlen2bootstrap();
703        get_rightson()->branchlen2bootstrap();
704    }
705}
706
707TreeNode *TreeNode::fixDeletedSon() {
708    // fix node after one son has been deleted
709    TreeNode *result = NULL;
710
711    if (leftson) {
712        gb_assert(!rightson);
713        result  = get_leftson();
714        leftson = NULL;
715    }
716    else {
717        gb_assert(!leftson);
718        gb_assert(rightson);
719
720        result   = get_rightson();
721        rightson = NULL;
722    }
723
724    // now 'result' contains the lasting tree
725    result->father = father;
726
727    if (remark_branch && !result->remark_branch) { // rescue remarks if possible
728        result->remark_branch    = remark_branch;
729        remark_branch = NULL;
730    }
731    if (gb_node && !result->gb_node) { // rescue group if possible
732        result->gb_node = gb_node;
733        gb_node         = NULL;
734    }
735
736    if (!result->father) {
737        get_tree_root()->change_root(this, result);
738    }
739
740    gb_assert(!is_root_node());
741
742    forget_origin();
743    forget_relatives();
744    delete this;
745
746    return result;
747}
748
749const TreeNode *TreeNode::ancestor_common_with(const TreeNode *other) const {
750    if (this == other) return this;
751    if (is_anchestor_of(other)) return this;
752    if (other->is_anchestor_of(this)) return other;
753    return get_father()->ancestor_common_with(other->get_father());
754}
755
756// --------------------------------------------------------------------------------
757
758#ifdef UNIT_TESTS
759#ifndef TEST_UNIT_H
760#include <test_unit.h>
761#endif
762
763void TEST_tree_iterator() {
764    GB_shell  shell;
765    GBDATA   *gb_main = GB_open("TEST_trees.arb", "r");
766    {
767        GB_transaction  ta(gb_main);
768        TreeNode       *tree = GBT_read_tree(gb_main, "tree_removal", new SimpleRoot);
769
770        int leafs = GBT_count_leafs(tree);
771        TEST_EXPECT_EQUAL(leafs, 17);
772        TEST_EXPECT_EQUAL(leafs_2_edges(leafs, UNROOTED), 31);
773
774        int iter_steps = ARB_edge::iteration_count(leafs);
775        TEST_EXPECT_EQUAL(iter_steps, 62);
776
777        const ARB_edge start = rootEdge(tree->get_tree_root());
778
779        // iterate forward + count (until same edge reached)
780        int      count       = 0;
781        int      count_leafs = 0;
782        ARB_edge edge        = start;
783        do {
784            ARB_edge next = edge.next();
785            TEST_EXPECT(next.previous() == edge); // test reverse operation
786            edge = next;
787            ++count;
788            if (edge.is_edge_to_leaf()) ++count_leafs;
789        }
790        while (edge != start);
791        TEST_EXPECT_EQUAL(count, iter_steps);
792        TEST_EXPECT_EQUAL(count_leafs, leafs);
793
794        // iterate backward + count (until same edge reached)
795        count       = 0;
796        count_leafs = 0;
797        edge        = start;
798        do {
799            ARB_edge next = edge.previous();
800            TEST_EXPECT(next.next() == edge); // test reverse operation
801            edge = next;
802            ++count;
803            if (edge.is_edge_to_leaf()) ++count_leafs;
804        }
805        while (edge != start);
806        TEST_EXPECT_EQUAL(count, iter_steps);
807        TEST_EXPECT_EQUAL(count_leafs, leafs);
808
809        if (tree) {
810            gb_assert(tree->is_root_node());
811            destroy(tree);
812        }
813    }
814    GB_close(gb_main);
815}
816
817TEST_PUBLISH(TEST_tree_iterator);
818
819#endif // UNIT_TESTS
820
821// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.