TreeNode.cxx

Go to the documentation of this file.

// ================================================================ //
//                                                                  //
//   File      : TreeNode.cxx                                       //
//   Purpose   :                                                    //
//                                                                  //
//   Coded by Ralf Westram (coder@reallysoft.de) in December 2013   //
//   Institute of Microbiology (Technical University Munich)        //
//   http://www.arb-home.de/                                        //
//                                                                  //
// ================================================================ //

#include "TreeNode.h"
#include <arb_progress.h>
#include <arb_str.h>
#include <map>
#include <set>
#include <cmath> // needed with 4.4.3 (but not with 4.7.3)

// ------------------
//      TreeRoot

TreeRoot::~TreeRoot() {
    deleteWithNodes = false; // avoid recursive call of ~TreeRoot (obsolete?)
    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
    // Note: destroying nodes from here is impossible (leads to pure virtual call, as derived class instance of 'this' is already under destruction)
}

void TreeRoot::change_root(TreeNode *oldroot, TreeNode *newroot) {
    // Note: when oldroot == NULp -> better use announce_tree_constructed()

    rt_assert(rootNode == oldroot);
    rt_assert(implicated(newroot, !newroot->father));
    rootNode = newroot;

    if (oldroot && oldroot->get_tree_root() && !oldroot->is_inside(newroot)) oldroot->set_tree_root(NULp); // unlink from this
    if (newroot && newroot->get_tree_root() != this) newroot->set_tree_root(this); // link to this
}

// --------------------
//      TreeNode

#if defined(PROVIDE_TREE_STRUCTURE_TESTS)

Validity TreeNode::is_valid() const {
    rt_assert(knownNonNull(this));
    Validity valid;

    TreeRoot *troot = get_tree_root();
    if (troot) {
        if (is_leaf()) {
            valid = Validity(!rightson && !leftson, "leaf has son");
        }
        else {
            valid            = Validity(rightson && leftson, "inner node lacks sons");
            if (valid) valid = get_rightson()->is_valid();
            if (valid) valid = get_leftson()->is_valid();
        }
        if (father) {
            if (valid) valid = Validity(is_inside(get_father()), "node not inside father subtree");
            if (valid) valid = Validity(troot->get_root_node()->is_ancestor_of(this), "root is not nodes ancestor");
            if (valid) valid = Validity(get_father()->get_tree_root() == troot, "node and father have different TreeRoot");
        }
        else {
            if (valid) valid = Validity(troot->get_root_node() == this, "node has no father, but isn't root-node");
            if (valid) valid = Validity(!is_leaf(), "root-node is leaf"); // leaf@root (tree has to have at least 2 leafs)
            if (valid) valid = Validity(has_valid_root_remarks(), "root-node has invalid remarks");
        }
    }
    else { // removed node (may be incomplete)
        if (is_leaf()) {
            valid = Validity(!rightson && !leftson, "(removed) leaf has son");
        }
        else {
            if (rightson)         valid = get_rightson()->is_valid();
            if (leftson && valid) valid = get_leftson()->is_valid();
        }
        if (father) {
            if (valid) valid = Validity(is_inside(get_father()), "(removed) node not inside father subtree");
            if (valid) valid = Validity(get_father()->get_tree_root() == troot, "(removed) node and father have different TreeRoot");
        }
    }
    return valid;
}
#endif // PROVIDE_TREE_STRUCTURE_TESTS

void TreeNode::set_tree_root(TreeRoot *new_root) {
    if (tree_root != new_root) {
        tree_root = new_root;
        if (leftson) get_leftson()->set_tree_root(new_root);
        if (rightson) get_rightson()->set_tree_root(new_root);
    }
}

void TreeNode::reorder_subtree(TreeOrder mode) {
    static const char *smallest_leafname; // has to be set to the alphabetically smallest name (when function exits)

    if (is_leaf()) {
        smallest_leafname = name;
    }
    else {
        int leftsize  = get_leftson() ->get_leaf_count();
        int rightsize = get_rightson()->get_leaf_count();

        {
            bool big_at_top    = leftsize>rightsize;
            bool big_at_bottom = leftsize<rightsize;
            bool swap_branches = (mode&ORDER_BIG_DOWN) ? big_at_top : big_at_bottom;
            if (swap_branches) swap_sons();
        }

        TreeOrder lmode, rmode;
        if (mode & (ORDER_BIG_TO_EDGE|ORDER_BIG_TO_CENTER)) { // symmetric
            if (mode & ORDER_ALTERNATING) mode = TreeOrder(mode^(ORDER_BIG_TO_EDGE|ORDER_BIG_TO_CENTER));

            if (mode & ORDER_BIG_TO_CENTER) {
                lmode = TreeOrder(mode | ORDER_BIG_DOWN);
                rmode = TreeOrder(mode & ~ORDER_BIG_DOWN);
            }
            else {
                lmode = TreeOrder(mode & ~ORDER_BIG_DOWN);
                rmode = TreeOrder(mode | ORDER_BIG_DOWN);
            }
        }
        else { // asymmetric
            if (mode & ORDER_ALTERNATING) mode = TreeOrder(mode^ORDER_BIG_DOWN);

            lmode = mode;
            rmode = mode;
        }

        get_leftson()->reorder_subtree(lmode);
        const char *leftleafname = smallest_leafname;

        get_rightson()->reorder_subtree(rmode);
        const char *rightleafname = smallest_leafname;

        if (leftleafname && rightleafname) {
            int name_cmp = strcmp(leftleafname, rightleafname);
            if (name_cmp <= 0) {
                smallest_leafname = leftleafname;
            }
            else {
                smallest_leafname = rightleafname;
                if (leftsize == rightsize) { // if sizes of subtrees are equal and rightleafname<leftleafname -> swap branches
                    const char *smallest_leafname_save = smallest_leafname;

                    swap_sons();
                    get_leftson ()->reorder_subtree(lmode); rt_assert(strcmp(smallest_leafname, rightleafname)== 0);
                    get_rightson()->reorder_subtree(rmode); rt_assert(strcmp(smallest_leafname, leftleafname) == 0);

                    smallest_leafname = smallest_leafname_save;
                }
            }
        }
    }
    rt_assert(smallest_leafname);
}

void TreeNode::reorder_tree(TreeOrder mode) {
    compute_tree();
    reorder_subtree(mode);
}

void TreeNode::rotate_subtree() {
    if (!is_leaf()) {
        swap_sons();
        get_leftson()->rotate_subtree();
        get_rightson()->rotate_subtree();
    }
}

void TreeNode::keelOver(TreeNode *prev, TreeNode *next, double len) {
    if (leftson == prev) {
        leftson = next;
        leftlen = len;

        if (keeledOver) {
            if (inverseLeft) keeledOver = false;
        }
        else {
            keeledOver  = true;
            inverseLeft = true;
        }
    }
    else {
        rightson = next;
        rightlen = len;

        if (keeledOver) {
            if (!inverseLeft) keeledOver = false;
        }
        else {
            keeledOver  = true;
            inverseLeft = false;
        }
    }
    father = prev;
}

void TreeNode::set_root() {
    if (at_root()) return; // already root

    TreeNode *old_root    = get_root_node();
    TreeNode *old_brother = is_inside(old_root->get_leftson()) ? old_root->get_rightson() : old_root->get_leftson();

    rt_assert(old_root->has_valid_root_remarks());

    // move remark branches to top
    {
        // Note: no remark is lost here (duplicate removed from old root; new duplicate created at new root)
        SmartCharPtr remarkPtr;
        if (!is_leaf()) remarkPtr = get_remark_ptr();
        for (TreeNode *node = this; node->father; node = node->get_father()) {
            std::swap(node->remark_branch, remarkPtr);
        }
    }

    GBT_LEN old_root_len = old_root->leftlen + old_root->rightlen;

    // new node & this init
    old_root->leftson  = this;
    old_root->rightson = father; // will be set later

    if (father->leftson == this) {
        old_root->leftlen = old_root->rightlen = father->leftlen*.5;
    }
    else {
        old_root->leftlen = old_root->rightlen = father->rightlen*.5;
    }

    TreeNode *next = get_father()->get_father();
    TreeNode *prev = old_root;
    TreeNode *pntr = get_father();

    if (father->leftson == this)    father->leftson = old_root; // to set the flag correctly

    // loop from father to son of root, rotate tree
    while  (next->father) {
        double len = (next->leftson == pntr) ? next->leftlen : next->rightlen;

        pntr->keelOver(prev, next, len);

        prev = pntr;
        pntr = next;
        next = next->get_father();
    }
    // now 'next' points to the old root, which has been destroyed above
    // 'pntr' points to one former son-of-root (the one nearer to new root)
    //
    // pointer at oldroot
    // pntr == brother before old root == next

    pntr->keelOver(prev, old_brother, old_root_len);

    old_brother->father = pntr;
    father              = old_root;

    rt_assert(get_root_node() == old_root); // the root node itself remains unchanged (its sons change)
    rt_assert(old_root->has_valid_root_remarks());
}

TreeNode *TreeNode::findLeafNamed(const char *wantedName) {
    TreeNode *found = NULp;
    if (is_leaf()) {
        if (name && strcmp(name, wantedName) == 0) found = this;
    }
    else {
        found             = get_leftson()->findLeafNamed(wantedName);
        if (!found) found = get_rightson()->findLeafNamed(wantedName);
    }
    return found;
}

// ----------------------------
//      find_innermost_edge

class NodeLeafDistance {
    GBT_LEN downdist, updist;
    enum { NLD_NODIST = 0, NLD_DOWNDIST, NLD_BOTHDIST } state;

public:

    NodeLeafDistance()
        : downdist(-1.0),
          updist(-1.0),
          state(NLD_NODIST)
    {}

    GBT_LEN get_downdist() const { rt_assert(state >= NLD_DOWNDIST); return downdist; }
    void set_downdist(GBT_LEN DownDist) {
        if (state < NLD_DOWNDIST) state = NLD_DOWNDIST;
        downdist = DownDist;
    }

    GBT_LEN get_updist() const { rt_assert(state >= NLD_BOTHDIST); return updist; }
    void set_updist(GBT_LEN UpDist) {
        if (state < NLD_BOTHDIST) state = NLD_BOTHDIST;
        updist = UpDist;
    }

};

class EdgeFinder {
    std::map<TreeNode*, NodeLeafDistance> data; // maximum distance to farthest leaf

    ARB_edge innermost;
    double   min_distdiff; // abs diff between up- and downdiff

    GBT_LEN calc_distdiff(GBT_LEN d1, GBT_LEN d2) { return fabs(d1-d2); }

    void insert_tree(TreeNode *node) {
        if (node->is_leaf()) {
            data[node].set_downdist(0.0);
        }
        else {
            insert_tree(node->get_leftson());
            insert_tree(node->get_rightson());

            data[node].set_downdist(std::max(data[node->get_leftson()].get_downdist()+node->leftlen,
                                             data[node->get_rightson()].get_downdist()+node->rightlen));
        }
    }

    void findBetterEdge_sub(TreeNode *node) {
        TreeNode *father  = node->get_father();
        TreeNode *brother = node->get_brother();

        GBT_LEN len      = node->get_branchlength();
        GBT_LEN brothLen = brother->get_branchlength();

        GBT_LEN upDist   = std::max(data[father].get_updist(), data[brother].get_downdist()+brothLen);
        GBT_LEN downDist = data[node].get_downdist();

        {
            GBT_LEN edge_dd = calc_distdiff(upDist, downDist);
            if (edge_dd<min_distdiff) { // found better edge
                innermost    = ARB_edge(node, father);
                min_distdiff = edge_dd;
            }
        }

        data[node].set_updist(upDist+len);

        if (!node->is_leaf()) {
            findBetterEdge_sub(node->get_leftson());
            findBetterEdge_sub(node->get_rightson());
        }
    }

    void findBetterEdge(TreeNode *node) {
        if (!node->is_leaf()) {
            findBetterEdge_sub(node->get_leftson());
            findBetterEdge_sub(node->get_rightson());
        }
    }

public:
    EdgeFinder(TreeNode *rootNode)
        : innermost(rootNode->get_leftson(), rootNode->get_rightson()) // root-edge
    {
        insert_tree(rootNode);

        TreeNode *lson = rootNode->get_leftson();
        TreeNode *rson = rootNode->get_rightson();

        GBT_LEN rootEdgeLen = rootNode->leftlen + rootNode->rightlen;

        GBT_LEN lddist = data[lson].get_downdist();
        GBT_LEN rddist = data[rson].get_downdist();

        data[lson].set_updist(rddist+rootEdgeLen);
        data[rson].set_updist(lddist+rootEdgeLen);

        min_distdiff = calc_distdiff(lddist, rddist);

        findBetterEdge(lson);
        findBetterEdge(rson);
    }

    const ARB_edge& innermost_edge() const { return innermost; }
};

ARB_edge TreeRoot::find_innermost_edge() {
    EdgeFinder edgeFinder(get_root_node());
    return edgeFinder.innermost_edge();
}

// ------------------------
//      multifurcation

class TreeNode::LengthCollector {
    typedef std::map<TreeNode*,GBT_LEN> LengthMap;
    typedef std::set<TreeNode*>         NodeSet;

    LengthMap eliminatedParentLength;
    LengthMap addedParentLength;

public:
    void eliminate_parent_edge(TreeNode *node) {
        rt_assert(!node->is_root_node());
        eliminatedParentLength[node] += parentEdge(node).eliminate();
    }

    void add_parent_length(TreeNode *node, GBT_LEN addLen) {
        rt_assert(!node->is_root_node());
        addedParentLength[node] += addLen;
    }

    void independent_distribution(bool useProgress) {
        // step 2: (see caller)
        arb_progress *progress    = NULp;
        int           redistCount = 0;
        if (useProgress) progress = new arb_progress("Distributing eliminated lengths", eliminatedParentLength.size());

        while (!eliminatedParentLength.empty()) { // got eliminated lengths which need to be distributed
            for (LengthMap::iterator from = eliminatedParentLength.begin(); from != eliminatedParentLength.end(); ++from) {
                ARB_edge elimEdge = parentEdge(from->first);
                GBT_LEN  elimLen  = from->second;

                elimEdge.virtually_distribute_length(elimLen, *this);
                if (progress) ++*progress;
            }
            eliminatedParentLength.clear(); // all distributed!

            // handle special cases where distributed length is negative and results in negative destination branches.
            // Avoid generating negative dest. branchlengths by
            // - eliminating the dest. branch
            // - redistributing the additional (negative) length (may cause additional negative lengths on other dest. branches)

            NodeSet handled;
            for (LengthMap::iterator to = addedParentLength.begin(); to != addedParentLength.end(); ++to) {
                ARB_edge affectedEdge     = parentEdge(to->first);
                GBT_LEN  additionalLen    = to->second;
                double   effective_length = affectedEdge.length() + additionalLen;

                if (effective_length<=0.0) { // negative or zero
                    affectedEdge.set_length(effective_length);
                    eliminate_parent_edge(to->first); // adds entry to eliminatedParentLength and causes another additional loop
                    handled.insert(to->first);
                }
            }

            if (progress && !eliminatedParentLength.empty()) {
                delete progress;
                ++redistCount;
                progress = new arb_progress(GBS_global_string("Redistributing negative lengths (#%i)", redistCount), eliminatedParentLength.size());
            }

            // remove all redistributed nodes
            for (NodeSet::iterator del = handled.begin(); del != handled.end(); ++del) {
                addedParentLength.erase(*del);
            }
        }

        // step 3:
        for (LengthMap::iterator to = addedParentLength.begin(); to != addedParentLength.end(); ++to) {
            ARB_edge affectedEdge     = parentEdge(to->first);
            GBT_LEN  additionalLen    = to->second;
            double   effective_length = affectedEdge.length() + additionalLen;

            affectedEdge.set_length(effective_length);
        }

        if (progress) delete progress;
    }
};

GBT_LEN ARB_edge::adjacent_distance() const {

    if (is_edge_to_leaf()) return 0.0;
    return next().length_or_adjacent_distance() + counter_next().length_or_adjacent_distance();
}

void ARB_edge::virtually_add_or_distribute_length_forward(GBT_LEN len, TreeNode::LengthCollector& collect) const {
    rt_assert(!is_nan_or_inf(len));
    if (length() > 0.0) {
        collect.add_parent_length(son(), len);
    }
    else {
        if (len != 0.0) virtually_distribute_length_forward(len, collect);
    }
}


void ARB_edge::virtually_distribute_length_forward(GBT_LEN len, TreeNode::LengthCollector& collect) const {
    rt_assert(is_normal(len));
    rt_assert(!is_edge_to_leaf()); // cannot forward anything (nothing beyond leafs)

    ARB_edge e1 = next();
    ARB_edge e2 = counter_next();

    GBT_LEN d1 = e1.length_or_adjacent_distance();
    GBT_LEN d2 = e2.length_or_adjacent_distance();

    len = len/(d1+d2);

    e1.virtually_add_or_distribute_length_forward(len*d1, collect);
    e2.virtually_add_or_distribute_length_forward(len*d2, collect);
}

void ARB_edge::virtually_distribute_length(GBT_LEN len, TreeNode::LengthCollector& collect) const {
    ARB_edge backEdge = inverse();
    GBT_LEN len_fwd, len_bwd;
    {
        GBT_LEN dist_fwd = adjacent_distance();
        GBT_LEN dist_bwd = backEdge.adjacent_distance();
        GBT_LEN lenW     = len/(dist_fwd+dist_bwd);
        len_fwd          = lenW*dist_fwd;
        len_bwd          = lenW*dist_bwd;

    }

    if (is_normal(len_fwd)) virtually_distribute_length_forward(len_fwd, collect);
    if (is_normal(len_bwd)) backEdge.virtually_distribute_length_forward(len_bwd, collect);
}

void TreeNode::eliminate_and_collect(const multifurc_limits& below, LengthCollector& collect) {
    rt_assert(!is_root_node());
    if (!is_leaf() || below.applyAtLeafs) {
        double value;
        if (is_leaf()) {
            value = 100.0;
            goto hack; // @@@ remove applyAtLeafs from multifurc_limits (does that really make sense? rethink!)
        }
        switch (parse_bootstrap(value, GB_warning)) {
            case REMARK_NONE:
                value = 100.0;
                // fall-through
            case REMARK_BOOTSTRAP:
            hack:
                if (value<below.bootstrap && get_branchlength_unrooted()<below.branchlength) {
                    collect.eliminate_parent_edge(this);
                }
                break;

            case REMARK_OTHER: break;
        }
    }

    if (!is_leaf()) {
        get_leftson() ->eliminate_and_collect(below, collect);
        get_rightson()->eliminate_and_collect(below, collect);
    }
}

void ARB_edge::multifurcate() {
    TreeNode::LengthCollector collector;
    collector.eliminate_parent_edge(son());
    collector.independent_distribution(false);
}
void TreeNode::multifurcate() {
    rt_assert(father); // cannot multifurcate at root; call with son of root to multifurcate root-edge
    if (father) parentEdge(this).multifurcate();
}

void TreeNode::set_branchlength_preserving(GBT_LEN new_len) {
    GBT_LEN old_len = get_branchlength_unrooted();
    GBT_LEN change  = new_len-old_len;

    char *old_remark = is_leaf() ? NULp : nulldup(get_remark());

    // distribute the negative 'change' to neighbours:
    set_branchlength_unrooted(-change);
    multifurcate();

    set_branchlength_unrooted(new_len);
    if (old_remark) {
        use_as_remark(old_remark); // restore remark (was removed by multifurcate())
    }
#if defined(ASSERTION_USED)
    else {
        rt_assert(has_no_remark());
    }
#endif
}

void TreeNode::multifurcate_whole_tree(const multifurc_limits& below) {
    TreeNode        *root = get_root_node();
    LengthCollector  collector;
    arb_progress progress("Multifurcating tree");

    // step 1:
    progress.subtitle("Collecting branches to eliminate");
    root->get_leftson()->eliminate_and_collect(below, collector);
    root->get_rightson()->eliminate_and_collect(below, collector);
    // root-edge is handled twice by the above calls.
    // Unproblematic: first call will eliminate root-edge, so second call will do nothing.

    // step 2 and 3:
    collector.independent_distribution(true);
}

void TreeNode::remove_bootstrap() {
    remark_branch.setNull();
    if (!is_leaf()) {
        get_leftson()->remove_bootstrap();
        get_rightson()->remove_bootstrap();
    }
}

GB_ERROR TreeNode::apply_aci_to_remarks(const char *aci, const GBL_call_env& callEnv) {
    GB_ERROR error = NULp;
    if (!is_leaf()) {
        {
            char *new_remark = GB_command_interpreter_in_env(remark_branch.isSet() ? remark_branch.content() : "", aci, callEnv);
            if (!new_remark) {
                error = GB_await_error();
            }
            else {
                freeset(new_remark, GBS_trim(new_remark));
                if (!new_remark[0]) { // skip empty results
                    free(new_remark);
                    remark_branch.setNull();
                }
                else {
                    remark_branch = new_remark;
                }
            }
        }

        if (!error) error = get_leftson()->apply_aci_to_remarks(aci, callEnv);
        if (!error) error = get_rightson()->apply_aci_to_remarks(aci, callEnv);
    }

    return error;
}
void TreeNode::reset_branchlengths() {
    if (!is_leaf()) {
        leftlen = rightlen = DEFAULT_BRANCH_LENGTH;

        get_leftson()->reset_branchlengths();
        get_rightson()->reset_branchlengths();
    }
}

void TreeNode::scale_branchlengths(double factor) {
    if (!is_leaf()) {
        leftlen  *= factor;
        rightlen *= factor;

        get_leftson()->scale_branchlengths(factor);
        get_rightson()->scale_branchlengths(factor);
    }
}

GBT_LEN TreeNode::sum_child_lengths() const {
    if (is_leaf()) return 0.0;
    return
        leftlen +
        rightlen +
        get_leftson()->sum_child_lengths() +
        get_rightson()->sum_child_lengths();
}

void TreeNode::bootstrap2branchlen() {
    if (is_leaf()) {
        set_branchlength_unrooted(DEFAULT_BRANCH_LENGTH);
    }
    else {
        if (!is_root_node()) {
            double bootstrap;
            if (parse_bootstrap(bootstrap, GB_warning) != REMARK_BOOTSTRAP) {
                bootstrap = 100.0; // no bootstrap means "100%"
            }

            double len = bootstrap/100.0;
            set_branchlength_unrooted(len);
        }
        get_leftson()->bootstrap2branchlen();
        get_rightson()->bootstrap2branchlen();
    }
}
void TreeNode::branchlen2bootstrap() {
    if (!is_leaf()) {
        remove_remark();
        if (!is_root_node()) {
            set_bootstrap(get_branchlength_unrooted()*100.0);
        }
        get_leftson()->branchlen2bootstrap();
        get_rightson()->branchlen2bootstrap();
    }
#if defined(ASSERTION_USED)
    else {
        rt_assert(has_no_remark());
    }
#endif
}
void TreeNode::branchlenXbootstrap() {
    if (!is_leaf()) {
        if (!is_root_node()) {
            if (!is_son_of_root() || is_leftson()) { // Note: condition depends on order of recursive calls below!
                double bootstrap;
                if (parse_bootstrap(bootstrap, GB_warning) != REMARK_BOOTSTRAP) {
                    bootstrap = 100.0; // no bootstrap means "100%"
                }
                double len = get_branchlength_unrooted();

                set_branchlength_unrooted(bootstrap/100.0);
                set_bootstrap(len*100.0);
            }
        }

        get_leftson()->branchlenXbootstrap();
        get_rightson()->branchlenXbootstrap();
    }
}

GBT_RemarkType TreeNode::parse_bootstrap(double& bootstrap, WarningConsumer report) {
    // If this function reports a warning to the caller, it will also repair the incorrect TreeNode(s).
    // Reporting a warning will cause the tree to be saved to DB in ARB_NTREE, if it has been loaded from DB.

    rt_assert(!is_leaf());      // only inner nodes may have bootstraps
    rt_assert(!is_root_node()); // root nodes have no bootstrap values (because they do not have a parent branch)

    GBT_RemarkType type = parse_remark(remark_branch.content(), bootstrap);
    if (is_son_of_root()) {
        TreeNode       *brother      = get_brother();
        double          brothersBootstrap;
        GBT_RemarkType  brothersType = parse_remark(brother->remark_branch.content(), brothersBootstrap);

        // Only inspect bootstrap(s) of root-edge, if at least one brother has type REMARK_BOOTSTRAP.
        if (type == REMARK_BOOTSTRAP) {
            if (brother->is_leaf()) {
                report(GBS_global_string("Ignored invalid bootstrap '%s' at root edge leading to a leaf", remark_branch.content()));
                remove_remark();
                type = REMARK_NONE;
            }
            else {
                switch (brothersType) {
                    case REMARK_BOOTSTRAP: {
                        double averageBootstrap = (bootstrap + brothersBootstrap) / 2.0;
                        double diff             = std::abs(bootstrap-brothersBootstrap);

                        SmartCharPtr myOrgRemark       = remark_branch;
                        SmartCharPtr brothersOrgRemark = brother->remark_branch;

                        set_bootstrap(averageBootstrap);
                        type = parse_remark(remark_branch.content(), bootstrap);

                        const double REMARKABLE_DIFF = 2.0;
#if defined(ASSERTION_USED)
                        const double EPSILON = 0.001 + 0.5; // "huge" EPSILON required, because set_bootstrap rounds value (causing 0.5 diff)
                        double diffToAverage = std::abs(bootstrap-averageBootstrap);
                        rt_assert(diffToAverage<EPSILON);
                        rt_assert(REMARKABLE_DIFF>EPSILON);
#endif

                        if (diff>REMARKABLE_DIFF) { // only warn for reasonable difference
                            report(GBS_global_string("Root-edge has conflicting support values '%s' and '%s' -> using average '%s'",
                                                     myOrgRemark.content(),
                                                     brothersOrgRemark.content(),
                                                     remark_branch.content()));
                        }
                        break;
                    }
                    case REMARK_OTHER:
                    case REMARK_NONE: {
                        // brother lacks bootstrap
                        const char *what    = GBS_global_string("Using comment '%s' on both sides of root-edge", remark_branch.content());
                        const char *warning = brothersType == REMARK_NONE ? what : GBS_global_string("%s (had '%s')", what, brother->remark_branch.content());

                        report(warning);
                        // -> use bootstrap of 'this' for both sons of root:
                        brother->use_as_remark(remark_branch);
                        break;
                    }
                }
            }
        }
        else if (brothersType == REMARK_BOOTSTRAP) {
            rt_assert(!brother->is_leaf());
            type = brother->parse_bootstrap(bootstrap, report);
        }
        // @@@ correct root remarks here. see code in TREE_load: ../SL/TREE_READ/TreeRead.cxx@CORR_ROOT_REM
    }
#if defined(DEBUG)
    else {
        // avoid using parse_bootstrap while constructing tree
        bool under_construction = father && !father->father; // otherwise would have entered the if-branch above
        rt_assert(!under_construction);    // fail if this tree is under construction
    }
#endif

    return type;
}
void TreeNode::set_bootstrap(double bootstrap_percent) {
    rt_assert(!is_leaf());      // only inner nodes may have bootstraps
    rt_assert(!is_root_node()); // root nodes have no bootstrap values (because they do not have a parent branch)

    int bootstrap_rounded = int(bootstrap_percent+0.5);
    if (bootstrap_rounded == 100) remove_remark();
    else                          use_as_remark(GBS_global_string_copy("%i%%", bootstrap_rounded));

    if (is_son_of_root()) {
        TreeNode *bro = get_brother();
        if (!bro->is_leaf()) bro->use_as_remark(remark_branch);
    }
}

TreeNode *TreeNode::fixDeletedSon() {
    // fix node after one son has been deleted
    TreeNode *result = NULp;

    if (leftson) {
        gb_assert(!rightson);
        result  = get_leftson();
        leftson = NULp;
    }
    else {
        gb_assert(!leftson);
        gb_assert(rightson);

        result   = get_rightson();
        rightson = NULp;
    }

    // now 'result' contains the lasting tree
    result->father = father;

    // rescue remarks if possible
    if (!is_leaf() && !result->is_leaf()) {
        if (get_remark() && !result->get_remark()) {
            result->use_as_remark(get_remark_ptr());
            remove_remark();
        }
    }

    if (gb_node && !result->gb_node) { // rescue group if possible
        result->gb_node = gb_node;
        gb_node         = NULp;
    }

    if (!result->father) {
        get_tree_root()->change_root(this, result);
    }

    gb_assert(!is_root_node());

    forget_origin();
    forget_relatives();
    delete this;

    return result;
}

const TreeNode *TreeNode::ancestor_common_with(const TreeNode *other) const {
    if (this == other) return this;
    if (is_ancestor_of(other)) return this;
    if (other->is_ancestor_of(this)) return other;
    return get_father()->ancestor_common_with(other->get_father());
}

int TreeNode::count_clades() const {
    if (is_leaf()) return is_clade();
    return get_leftson()->count_clades() + get_rightson()->count_clades() + is_clade();
}

#if defined(ASSERTION_USED) || defined(PROVIDE_TREE_STRUCTURE_TESTS)
bool TreeNode::has_valid_root_remarks() const {
    // tests whether the root-edge has a valid remark:
    // - if one son is a leaf, the other son may contain the remark for the root-edge
    // - if no son is a leaf, both sons shall contain the same remark (or none)
    rt_assert(is_root_node());
    return implicated(!get_leftson()->is_leaf() && !get_rightson()->is_leaf(),
                      ARB_strNULLcmp(get_leftson()->get_remark(), get_rightson()->get_remark()) == 0);
}

#endif

// --------------------------------------------------------------------------------

#ifdef UNIT_TESTS
#ifndef TEST_UNIT_H
#include <test_unit.h>
#endif

void TEST_tree_iterator() {
    GB_shell  shell;
    GBDATA   *gb_main = GB_open("TEST_trees.arb", "r");
    {
        GB_transaction  ta(gb_main);
        TreeNode       *tree = GBT_read_tree(gb_main, "tree_removal", new SimpleRoot);

        int leafs = GBT_count_leafs(tree);
        TEST_EXPECT_EQUAL(leafs, 17);
        TEST_EXPECT_EQUAL(leafs_2_edges(leafs, UNROOTED), 31);

        int iter_steps = ARB_edge::iteration_count(leafs);
        TEST_EXPECT_EQUAL(iter_steps, 62);

        // test edge-cases of iteration_count:
        {
            TEST_EXPECT_EQUAL(ARB_edge::iteration_count(3), 6); // 3 leafs -> 4 nodes -> 3 edges in 2 directions -> 6 iterations
            TEST_EXPECT_EQUAL(ARB_edge::iteration_count(2), 2); // 2 leafs -> 2 nodes -> 1 edge in 2 directions -> 2 iterations
            TEST_EXPECT_EQUAL(ARB_edge::iteration_count(1), 0); // one leaf -> invalid tree -> will not iterate -> return 0
            TEST_EXPECT_EQUAL(ARB_edge::iteration_count(0), 0); // no leafs -> invalid tree -> will not iterate -> return 0
        }

        const ARB_edge start = rootEdge(tree->get_tree_root());

        // iterate forward + count (until same edge reached)
        int      count       = 0;
        int      count_leafs = 0;
        ARB_edge edge        = start;
        do {
            ARB_edge next = edge.next();
            TEST_EXPECT(next.previous() == edge); // test reverse operation
            edge = next;
            ++count;
            if (edge.is_edge_to_leaf()) ++count_leafs;
        }
        while (edge != start);
        TEST_EXPECT_EQUAL(count, iter_steps);
        TEST_EXPECT_EQUAL(count_leafs, leafs);

        // iterate backward + count (until same edge reached)
        count       = 0;
        count_leafs = 0;
        edge        = start;
        do {
            ARB_edge next = edge.previous();
            TEST_EXPECT(next.next() == edge); // test reverse operation
            edge = next;
            ++count;
            if (edge.is_edge_to_leaf()) ++count_leafs;
        }
        while (edge != start);
        TEST_EXPECT_EQUAL(count, iter_steps);
        TEST_EXPECT_EQUAL(count_leafs, leafs);

        if (tree) {
            gb_assert(tree->is_root_node());
            destroy(tree);
        }
    }
    GB_close(gb_main);
}

void TEST_tree_branch_modifications() {
    GB_shell  shell;
    GBDATA   *gb_main = GB_open("TEST_trees.arb", "r");
    {
        const char *newick_bs2bl = "(((((((CloTyro3:.1,CloTyro4:.1):.4,CloTyro2:.1):0,CloTyrob:.1):.97,CloInnoc:.1):0,CloBifer:.1):.53,(((CloButy2:.1,CloButyr:.1):1,CloCarni:.1):.33,CloPaste:.1):.97):.5,((((CorAquat:.1,CurCitre:.1):1,CorGluta:.1):.17,CelBiazo:.1):.4,CytAquat:.1):.5);";                                                                // bootstrap2branchlen
        const char *newick_bl2bs = "(((((((CloTyro3,CloTyro4)'3%',CloTyro2)'2%',CloTyrob)'27%',CloInnoc)'6%',CloBifer)'12%',(((CloButy2,CloButyr)'56%',CloCarni)'1%',CloPaste)'13%')'16%',((((CorAquat,CurCitre)'10%',CorGluta)'5%',CelBiazo)'21%',CytAquat)'16%');";                                                                                                                                                                            // branchlen2bootstrap
        const char *newick       = "(((((((CloTyro3:1.046,CloTyro4:.061)'40%':.026,CloTyro2:.017)'0%':.017,CloTyrob:.009)'97%':.274,CloInnoc:.371)'0%':.057,CloBifer:.388)'53%':.124,(((CloButy2:.009,CloButyr:0):.564,CloCarni:.12)'33%':.01,CloPaste:.179)'97%':.131):.081,((((CorAquat:.084,CurCitre:.058):.103,CorGluta:.522)'17%':.053,CelBiazo:.059)'40%':.207,CytAquat:.711):.081);";

#define LEFT_SWAP_EVEN  "((((((CloTyro3:1.046,CloTyro4:.061)'40%':.03,CloTyro2:.017)'0%':.02,CloTyrob:.009)'97%':.27,CloInnoc:.371)'0%':.06,CloBifer:.388)'53%':.12,(((CloButy2:.009,CloButyr:0):.56,CloCarni:.12)'33%':.01,CloPaste:.179)'97%':.13)"
#define LEFT_SWAP_ODD   "((((((CloTyro3:1.046,CloTyro4:.061)'3%':.4,CloTyro2:.017)'2%':0,CloTyrob:.009)'27%':.97,CloInnoc:.371)'6%':0,CloBifer:.388)'12%':.53,(((CloButy2:.009,CloButyr:0)'56%':1,CloCarni:.12)'1%':.33,CloPaste:.179)'13%':.97)"
#define RIGHT_SWAP_EVEN "((((CorAquat:.084,CurCitre:.058):.1,CorGluta:.522)'17%':.05,CelBiazo:.059)'40%':.21,CytAquat:.711)"
#define RIGHT_SWAP_ODD  "((((CorAquat:.084,CurCitre:.058)'10%':1,CorGluta:.522)'5%':.17,CelBiazo:.059)'21%':.4,CytAquat:.711)"

        const char *newick_blXbs_2 = "(" LEFT_SWAP_EVEN ":.08," RIGHT_SWAP_EVEN ":.08);"; // branchlenXbootstrap
        const char *newick_blXbs_1 = "(" LEFT_SWAP_ODD   "'16%':.5," RIGHT_SWAP_ODD  "'16%':.5);"; // branchlenXbootstrap

        const char *right_newick            = "((((CorAquat:.084,CurCitre:.058):.103,CorGluta:.522)'17%':.053,CelBiazo:.059)'40%':.207,CytAquat:.711):.081;";
        const char *right_newick_preserved1 = "((((CorAquat:.084,CurCitre:.058):.103,CorGluta:.522)'17%':.06,CelBiazo:.029)'40%':.231,CytAquat:.711):.081;"; // set_branchlength_preserving to 0.029 at CelBiazo
        const char *right_newick_preserved2 = "((((CorAquat:.084,CurCitre:.058):.103,CorGluta:.522)'17%':.041,CelBiazo:.117)'40%':.161,CytAquat:.711):.081;"; // set_branchlength_preserving to 0.117 at CelBiazo

        NewickFormat BSLEN = NewickFormat(nREMARK|nLENGTH);

        GB_transaction ta(gb_main);

        // test bootstrap2branchlen:
        {
            TreeNode *tree = GBT_read_tree(gb_main, "tree_test", new SimpleRoot);
            TEST_EXPECT_NEWICK(BSLEN, tree, newick);
            gb_assert(tree->has_valid_root_remarks());
            tree->bootstrap2branchlen();
            TEST_EXPECT_NEWICK(nLENGTH, tree, newick_bs2bl);
            gb_assert(tree->has_valid_root_remarks());
            destroy(tree);
        }

        // test branchlen2bootstrap:
        {
            TreeNode *tree = GBT_read_tree(gb_main, "tree_test", new SimpleRoot);
            TEST_EXPECT_NEWICK(BSLEN, tree, newick);
            tree->branchlen2bootstrap();
            TEST_EXPECT_NEWICK(nREMARK, tree, newick_bl2bs);
            gb_assert(tree->has_valid_root_remarks());
            destroy(tree);
        }
        // test branchlenXbootstrap:
        {
            TreeNode *tree = GBT_read_tree(gb_main, "tree_test", new SimpleRoot);
            TEST_EXPECT_NEWICK(BSLEN, tree, newick);

            TEST_EXPECT_DIFFERENT(newick_blXbs_2, newick); // did change due to rounding issues

            for (int repeat = 1; repeat<=3; ++repeat) {
                // swap operation stabilizes after 2 swaps

                tree->branchlenXbootstrap(); // 1st swap
                TEST_EXPECT_NEWICK(BSLEN, tree, newick_blXbs_1);
                gb_assert(tree->has_valid_root_remarks());

                tree->branchlenXbootstrap(); // 2nd swap
                TEST_EXPECT_NEWICK(BSLEN, tree, newick_blXbs_2);
                gb_assert(tree->has_valid_root_remarks());
            }

            destroy(tree);
        }

        {
            TreeNode *tree = GBT_read_tree(gb_main, "tree_test", new SimpleRoot);

            TreeNode *right    = tree->get_rightson();
            TreeNode *CelBiazo = right->findLeafNamed("CelBiazo");

            gb_assert(tree->has_valid_root_remarks());

            TEST_EXPECT_NEWICK(BSLEN, right, right_newick);

            CelBiazo->set_branchlength_preserving(CelBiazo->get_branchlength() * 0.5);
            TEST_EXPECT_NEWICK(BSLEN, right, right_newick_preserved1);

            CelBiazo->set_branchlength_preserving(CelBiazo->get_branchlength() * 4.0);
            TEST_EXPECT_NEWICK(BSLEN, right, right_newick_preserved2);

            gb_assert(tree->has_valid_root_remarks());

            destroy(tree);
        }
    }
    GB_close(gb_main);
}

TEST_PUBLISH(TEST_tree_branch_modifications);

#endif // UNIT_TESTS

// --------------------------------------------------------------------------------