source: branches/stable/ARBDB/TreeNode.cxx

// ================================================================ //
//                                                                  //
//   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) {
    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) {
    /*! beautify tree (does not change topology, only swaps branches)
     */
    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) {
    /*! atomic part of set_root operation that keels over edges between new and old root
     * @param prev previously a son of 'this', will become father
     * @param next previously the father of 'this', will become son
     * @param len  length of branch to nextSon
     */

    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() {
    /*! set the root at parent edge of this
     * pointers to tree-nodes remain valid, but all parent-nodes of this change their meaning
     * (afterwards they will point to [father+brother] instead of [this+brother])
     * esp. pointers to the root-node will still point to the root-node (which only changed children)
     */

    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 {
    //! return length of edges starting from this->dest()

    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 {
    /*! distribute length to edges adjacent in edge direction (i.e. edges starting from this->dest())
     * Split 'len' proportional to adjacent edges lengths.
     *
     * Note: length will not be distributed to tree-struction itself (yet), but collected in 'collect'
     */

    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 {
    /*! distribute length to all adjacent edges.
     * Longer edges receive more than shorter ones.
     *
     * Edges with length zero will not be changed, instead both edges "beyond"
     * the edge will be affected (they will be affected equally to direct edges,
     * because edges at multifurcations are considered to BE direct edges).
     *
     * Note: length will not be distributed to tree-struction itself (yet), but collected in 'collect'
     */

    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) {
    /*! eliminate edges specified by 'below' and
     * store their length in 'collect' for later distribution.
     */
    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)) {
            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() {
    /*! eliminate edge and distribute length to adjacent edges
     * - sets its length to zero
     * - removes remark (bootstrap)
     * - distributes length to neighbour-branches
     */
    TreeNode::LengthCollector collector;
    collector.eliminate_parent_edge(son());
    collector.independent_distribution(false);
}
void TreeNode::multifurcate() {
    /*! eliminate branch from 'this' to 'father' (or brother @ root)
     * @see ARB_edge::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) {
    /*! set branchlength to 'new_len' while preserving overall distance in tree.
     *
     * Always works on unrooted tree (i.e. lengths @ root are treated correctly).
     * Length is preserved as in multifurcate()
     */

    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) {
    /*! multifurcate all branches specified by 'below'
     * - step 1: eliminate all branches, store eliminated lengths
     * - step 2: calculate length distribution for all adjacent branches (proportionally to original length of each branch)
     * - step 3: apply distributed length
     */
    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() {
    //! copy bootstraps to branchlengths
    if (is_leaf()) {
        set_branchlength_unrooted(DEFAULT_BRANCH_LENGTH);
    }
    else {
        if (father) {
            double         bootstrap;
            GBT_RemarkType rtype = parse_bootstrap(bootstrap);

            if (rtype == REMARK_BOOTSTRAP) {
                double len = bootstrap/100.0;
                set_branchlength_unrooted(len);
            }
            else {
                set_branchlength_unrooted(1.0); // no bootstrap means "100%"
            }
        }
        get_leftson()->bootstrap2branchlen();
        get_rightson()->bootstrap2branchlen();
    }
}

void TreeNode::branchlen2bootstrap() {
    //! copy branchlengths to bootstraps
    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
}

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);

        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 *left_newick       = "((((((CloTyro3:1.046,CloTyro4:0.061)'40%':0.026,CloTyro2:0.017)'0%':0.017,CloTyrob:0.009)'97%':0.274,CloInnoc:0.371)'0%':0.057,CloBifer:0.388)'53%':0.124,(((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120)'33%':0.010,CloPaste:0.179)'97%':0.131):0.081;";
        const char *left_newick_bs2bl = "((((((CloTyro3:0.100,CloTyro4:0.100):0.400,CloTyro2:0.100):0.000,CloTyrob:0.100):0.970,CloInnoc:0.100):0.000,CloBifer:0.100):0.530,(((CloButy2:0.100,CloButyr:0.100):1.000,CloCarni:0.100):0.330,CloPaste:0.100):0.970):0.500;"; // bootstrap2branchlen

        const char *right_newick            = "((((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522)'17%':0.053,CelBiazo:0.059)'40%':0.207,CytAquat:0.711):0.081;";
        const char *right_newick_preserved1 = "((((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522)'17%':0.060,CelBiazo:0.029)'40%':0.231,CytAquat:0.711):0.081;"; // set_branchlength_preserving to 0.029 at CelBiazo
        const char *right_newick_preserved2 = "((((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522)'17%':0.041,CelBiazo:0.117)'40%':0.161,CytAquat:0.711):0.081;"; // set_branchlength_preserving to 0.117 at CelBiazo
        const char *right_newick_bl2bs      = "((((CorAquat,CurCitre)'10%',CorGluta)'5%',CelBiazo)'21%',CytAquat)'100%';"; // branchlen2bootstrap

        NewickFormat BSLEN = NewickFormat(nREMARK|nLENGTH);

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

            TreeNode *left  = tree->get_leftson();
            TreeNode *right = tree->get_rightson();

            TEST_EXPECT_NEWICK(BSLEN, left, left_newick);
            TEST_EXPECT_NEWICK(BSLEN, right, right_newick);

            left->bootstrap2branchlen();
            right->branchlen2bootstrap();

            TEST_EXPECT_NEWICK(nLENGTH, left,  left_newick_bs2bl);
            TEST_EXPECT_NEWICK(nREMARK, right, right_newick_bl2bs);

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

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

            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);

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

TEST_PUBLISH(TEST_tree_branch_modifications);

#endif // UNIT_TESTS

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