placement/function/tree.cpp

Go to the documentation of this file.

/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2022 Lucas Czech
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
    Contact:
    Lucas Czech <lczech@carnegiescience.edu>
    Department of Plant Biology, Carnegie Institution For Science
    260 Panama Street, Stanford, CA 94305, USA
*/
 
#include "genesis/placement/function/tree.hpp"
 
#include "genesis/placement/function/helper.hpp"
#include "genesis/tree/common_tree/operators.hpp"
#include "genesis/tree/common_tree/tree.hpp"
#include "genesis/tree/function/functions.hpp"
#include "genesis/tree/function/manipulation.hpp"
#include "genesis/tree/function/operators.hpp"
#include "genesis/tree/tree.hpp"
 
#include <algorithm>
#include <cassert>
#include <numeric>
#include <stdexcept>
#include <utility>
#include <vector>
 
namespace genesis {
namespace placement {
 
// =================================================================================================
//     Local Helper Functions
// =================================================================================================
 
// Helper struct that stores a placement along with its containing pquery.
struct LabelledTreePlacementPair
{
    LabelledTreePlacementPair( Pquery const* pq, PqueryPlacement const* pp )
        : pquery(    pq )
        , placement( pp )
    {}
 
    Pquery const*          pquery;
    PqueryPlacement const* placement;
};
 
// -----------------------------------------------------
//     Placement Pairs per Edge
// -----------------------------------------------------
 
// For each edge, get a list of placements - but only the most probable ones per Pquery.
// This is similar to placements_per_edge(), but we also need a pointer to the pquery, in
// order to get its name(s).
std::vector< std::vector< LabelledTreePlacementPair >> labelled_tree_placement_pairs_per_edge_(
    Sample const& sample
) {
    std::vector< std::vector< LabelledTreePlacementPair >> place_map;
    place_map.resize( sample.tree().edge_count() );
    for( auto const& pqry : sample.pqueries() ) {
 
        // Find the most probable placement.
        auto max_it = std::max_element(
            pqry.placements().begin(),
            pqry.placements().end(),
            []( PqueryPlacement const& lhs, PqueryPlacement const& rhs ){
                return lhs.like_weight_ratio < rhs.like_weight_ratio;
            }
        );
 
        // If there is one (i.e., the pquery has at least one placement),
        // add it to the list for its edge.
        if( max_it != pqry.placements().end() ) {
            place_map[ max_it->edge().index() ].emplace_back( &pqry, &*max_it );
        }
 
        // Alternative handcrafted version...
        // PqueryPlacement const* max_p = nullptr;
        // double max_v = std::numeric_limits<double>::lowest();
        // for( auto const& place : pqry.placements() ) {
        //     if( place.like_weight_ratio > max_v ) {
        //         max_v = place.like_weight_ratio;
        //         max_p = &place;
        //     }
        // }
        // if( max_p ) {
        //     place_map[ max_p->edge().index() ].emplace_back( &pqry, max_p );
        // }
    }
    return place_map;
}
 
// -----------------------------------------------------
//     Single Placement on the Edge
// -----------------------------------------------------
 
// Helper function that adds one placement, if it is the only one on its edge.
void labelled_tree_add_lonely_placement_(
    tree::Tree&          tree,
    tree::TreeEdge&      edge,
    LabelledTreePlacementPair const& placement_pair,
    std::string const& name_prefix
) {
    // Add the new edges to the tree and get all necessary edges.
    auto& pendant_node  = tree::add_new_leaf_node( tree, edge );
    auto& pendant_edge  = pendant_node.link().edge();
    auto& proximal_edge = pendant_edge.primary_link().next().edge();
    auto& distal_edge   = pendant_edge.primary_link().next().next().edge();
 
    // The inner node is new, so it should be bifurcating. The other one is a leaf.
    assert( degree( pendant_edge.primary_node() ) == 3 );
    assert( is_leaf( pendant_edge.secondary_node() ));
 
    // Some shorthands to the edge data.
    auto& pend_data = pendant_edge.data<tree::CommonEdgeData>();
    auto& prox_data = proximal_edge.data<tree::CommonEdgeData>();
    auto& dist_data  = distal_edge.data<tree::CommonEdgeData>();
 
    // The pendant and distal edges are new, so they should have default branch lengths.
    assert( pend_data.branch_length == 0.0 );
    assert( dist_data.branch_length == 0.0 );
 
    // Set all three branch lengths. Make sure they are not negative.
    double const distal_len = prox_data.branch_length - placement_pair.placement->proximal_length;
    pend_data.branch_length = placement_pair.placement->pendant_length;
    prox_data.branch_length = placement_pair.placement->proximal_length;
    dist_data.branch_length = std::max( 0.0, distal_len );
 
    // Set the leaf node name, if there is one. Also, assert that ther is at most one -
    // otherwise, this is not a lonely placement any more!
    assert( placement_pair.pquery->name_size() <= 1 );
    if( placement_pair.pquery->name_size() == 1 ) {
        auto& pendant_node_data = pendant_node.data<tree::CommonNodeData>();
        pendant_node_data.name = name_prefix + placement_pair.pquery->name_at(0).name;
    }
}
 
// -----------------------------------------------------
//     Multiple Placements, Fully Resolved
// -----------------------------------------------------
 
// Helper function that adds all placements for one edge, fully resolved.
void labelled_tree_process_edge_fully_resolved_(
    tree::Tree&                       tree,
    tree::TreeEdge&                   edge,
    std::vector<LabelledTreePlacementPair> const& placement_pairs,
    std::string const& name_prefix
) {
    // In each step, we add a new node, along the original edge, splitting it more and more.
    // The branch lengths of those fragments are the differences between the proximal lengths
    // of the placements.
    // To handle this, store the proximal length that we already used along the branch.
    double used_length = 0.0;
 
    // Also, we need the original branch length, because it will be lost in the process,
    // but later needed.
    auto const orig_length = edge.data<tree::CommonEdgeData>().branch_length;
 
    // Store an edge pointer at which we insert the next placement edge.
    auto insertion_edge = &edge;
 
    // Process all placements.
    for( auto const& placement_pair : placement_pairs ) {
        // Each name gets its own branch.
        for( auto const& pquery_name : placement_pair.pquery->names() ) {
 
            // Create the new edges.
            auto& pendant_node  = tree::add_new_leaf_node( tree, *insertion_edge );
            auto& pendant_edge  = pendant_node.link().edge();
            auto& proximal_edge = pendant_edge.primary_link().next().edge();
            auto& distal_edge   = pendant_edge.primary_link().next().next().edge();
 
            // The primary node is new, so it should be bifurcating and a leaf.
            assert( degree( pendant_edge.primary_node() ) == 3 );
            assert( is_leaf( pendant_edge.secondary_node() ));
 
            // Some shorthands to the edge data. We dont need dist data; its branch length
            // will be set in a later step. Either it becomes the prox edge in the next
            // iteration, or it is the final edge after the loop.
            auto& pend_data = pendant_edge.data<tree::CommonEdgeData>();
            auto& prox_data = proximal_edge.data<tree::CommonEdgeData>();
 
            // The pendant and distal edges are new, so they should have default branch lengths.
            assert( pend_data.branch_length == 0.0 );
            assert( distal_edge.data<tree::CommonEdgeData>().branch_length == 0.0 );
 
            // We sorted the placements by prox len. Thus, this should hold.
            assert( placement_pair.placement->proximal_length >= used_length );
 
            // Set branch properties.
            pend_data.branch_length = placement_pair.placement->pendant_length;
            prox_data.branch_length = placement_pair.placement->proximal_length - used_length;
 
            // Set the leaf name.
            auto& node_data = pendant_node.data<tree::CommonNodeData>();
            node_data.name = name_prefix + pquery_name.name;
 
            // Set new values for the keeping-track variables, for the next iteration.
            used_length    = placement_pair.placement->proximal_length;
            insertion_edge = &distal_edge;
        }
    }
 
    // Now, the insertion edge is the last remaining edge part towards the original
    // end of the edge. Its branch length is not yet set. Assert this, then set it to what
    // is left of it.
    // If the remaining length is however smaller than the original one, that means we have
    // placements on this edge with a proximal length greater than the original edge length.
    // This is a bit weird in terms of placement, but might happen due to BLO in the placement.
    // In this case, we want to avoid negative length, so simply set it to zero.
    // That means that in total the edge now grew due to its placements. So the tree is not
    // optimized any more - but for the purposes of this function, this should be okay.
    auto& end_edge_data = insertion_edge->data<tree::CommonEdgeData>();
    assert( end_edge_data.branch_length == 0.0 );
    if( used_length < orig_length ) {
        end_edge_data.branch_length = orig_length - used_length;
    }
}
 
// -----------------------------------------------------
//     Multiple Placements, Multifurcating
// -----------------------------------------------------
 
// Helper function that adds all placements for one edge, multifurcating,
// if there is more than one placement for the edge.
void labelled_tree_process_edge_multifurcating_(
    tree::Tree&                       tree,
    tree::TreeEdge&                   edge,
    std::vector<LabelledTreePlacementPair> const& placement_pairs,
    std::string const& name_prefix
) {
    // Add a new leaf node to the tree, attached to the middle of the given edge (this spltis
    // the edge in half and adds one more node there). This will be the base to which we then
    // attach all further nodes, multifurcating.
    auto& new_node  = tree::add_new_leaf_node( tree, edge );
    auto& base_edge = new_node.link().edge();
    auto& pri_edge  = base_edge.primary_link().next().edge();
    auto& sec_edge  = base_edge.primary_link().next().next().edge();
 
    // The base node is new, so it should be bifurcating and a leaf
    assert( degree( base_edge.primary_node() ) == 3 );
    assert( is_leaf( base_edge.secondary_node() ));
 
    // Some shorthands to the edge data.
    auto& base_data = base_edge.data<tree::CommonEdgeData>();
    auto& pri_data  = pri_edge.data<tree::CommonEdgeData>();
    auto& sec_data  = sec_edge.data<tree::CommonEdgeData>();
 
    // The pendant and distal edges are new, so they should have default branch lengths.
    assert( base_data.branch_length == 0.0 );
    assert( sec_data.branch_length  == 0.0 );
 
    // Calculate the average attachment position (proximal length) of all placements.
    auto const avg_prox_len = std::accumulate(
        placement_pairs.begin(),
        placement_pairs.end(),
        0.0,
        []( double in, LabelledTreePlacementPair const& cur ){
            return in + cur.placement->proximal_length;
        }
    ) / static_cast<double>( placement_pairs.size() );
 
    // Set the branch lengths of the two splitted edges (the ones that were one edge before)
    // so that they each take their share of the original branch length, divided using
    // the average proximal length of the placements. This way, the multifurcation that contains
    // the placements will sit at the sweet spot.
    // Set the second one first, as it uses the former length of the first one (i.e., the
    // original branch length).
    // Also, make sure that each new length is non-negative. This might happen if
    // the placements are somehow weirdly placed.
    sec_data.branch_length = std::max( 0.0, pri_data.branch_length - avg_prox_len );
    pri_data.branch_length = std::max( 0.0, avg_prox_len );
 
    // This function should only be called with at least one placement. Assert this.
    // Then, we can safely get the placement with the smallest pendant length.
    assert( placement_pairs.size() > 0 );
    auto const min_pen_len = std::max( 0.0, std::min_element(
        placement_pairs.begin(),
        placement_pairs.end(),
        []( LabelledTreePlacementPair const& lhs, LabelledTreePlacementPair const& rhs ){
            return lhs.placement->pendant_length < rhs.placement->pendant_length;
        }
    )->placement->pendant_length );
 
    // TODO we sorted them first, so taking the first placement should also work:
    // (if so, we can drop the above search)
    assert( min_pen_len == placement_pairs[0].placement->pendant_length );
 
    // Use the min bl for the multifurcation base edge.
    base_data.branch_length = min_pen_len;
 
    // Now add all placements as edges to the base node (secondary node of the base edge).
    for( auto const& placement_pair : placement_pairs ) {
        // Each name gets its own branch.
        for( auto const& pquery_name : placement_pair.pquery->names() ) {
 
            // Make a new leaf node for this placement.
            auto& p_node = tree::add_new_node( tree, new_node );
            auto& p_edge = p_node.link().edge();
 
            // Set the pendant branch length. It is subtracted by
            // the min pen length, as this is already incorporated into the base's bl.
            auto& p_data = p_edge.data<tree::CommonEdgeData>();
            assert( placement_pair.placement->pendant_length >= min_pen_len );
            p_data.branch_length = placement_pair.placement->pendant_length - min_pen_len;
 
            // Set the leaf node name.
            auto& p_node_data = p_node.data<tree::CommonNodeData>();
            p_node_data.name = name_prefix + pquery_name.name;
        }
    }
}
 
// =================================================================================================
//     Helper Functions
// =================================================================================================
 
tree::Tree labelled_tree(
    Sample const&      sample,
    bool const         fully_resolve,
    std::string const& name_prefix
) {
    // Get a copy of the original tree that contains only common data.
    return labelled_tree(
        sample,
        tree::convert_to_common_tree( sample.tree() ),
        fully_resolve,
        name_prefix
    );
}
 
tree::Tree labelled_tree(
    Sample const&      sample,
    tree::Tree const&  tree,
    bool               fully_resolve,
    std::string const& name_prefix
) {
    // -----------------------------------------------------
    //     Preparation
    // -----------------------------------------------------
 
    // Get a copy of the original tree that we can add edges to.
    auto result = tree;
 
    // Check whether the tree is compatible with the sample tree.
    // This is a bit wasteful for cases when this function is called directly from its sample
    // version above (then, we can be sure that it is compatible), but necessary otherwise.
    if( ! tree::identical_topology( result, sample.tree() )) {
        throw std::runtime_error(
            "Tree provided for labelled_tree() is not compatible "
            "with the Tree of the provided Sample."
        );
    }
 
    // Get a list of the edge pointers so that we can iterate and add to the tree at the same time.
    std::vector< tree::TreeEdge* > edge_list;
    edge_list.reserve( result.edge_count() );
    for( auto& edge : result.edges() ) {
        edge_list.push_back( &edge );
    }
    assert( edge_list.size() == sample.tree().edge_count() );
 
    auto place_map = labelled_tree_placement_pairs_per_edge_( sample );
    assert( place_map.size() == edge_list.size() );
 
    // Sort the placements per edge, so that we can attach them to the tree in order.
    // For a fully resolved tree, we want to sort by proximal length, as this is the order in which
    // the placements will form new nodes on the tree.
    // For a multifurcating tree, we sort by pendant length, so that the short ones come first.
    auto sort_placements = fully_resolve
        ? []( LabelledTreePlacementPair const& lhs, LabelledTreePlacementPair const& rhs ){
              return lhs.placement->proximal_length < rhs.placement->proximal_length;
          }
        : []( LabelledTreePlacementPair const& lhs, LabelledTreePlacementPair const& rhs ){
              return lhs.placement->pendant_length < rhs.placement->pendant_length;
          }
    ;
    for( auto& edge_l : place_map ) {
        std::sort( edge_l.begin(), edge_l.end(), sort_placements );
    }
 
    // -----------------------------------------------------
    //     Main Loop over Edges
    // -----------------------------------------------------
 
    // Process each edge.
    for( size_t edge_i = 0; edge_i < edge_list.size(); ++edge_i ) {
        // Some safety. There needs to be a valid pointer to an edge.
        assert( edge_list[ edge_i ] );
 
        // Nothing to do if there are no placements for this edge.
        if( place_map[ edge_i ].size() == 0 ) {
            continue;
        }
 
        // If there is only one placement with at most one name,
        // both algorithms behave the same, so shortcut this.
        if( place_map[ edge_i ].size() == 1 && place_map[ edge_i ][0].pquery->name_size() <= 1 ) {
            assert( place_map[ edge_i ][0].pquery && place_map[ edge_i ][0].placement );
            labelled_tree_add_lonely_placement_(
                result, *edge_list[ edge_i ], place_map[ edge_i ][0], name_prefix
            );
            continue;
        }
 
        // Select which algorithm to use for more than one placement at the edge.
        if( fully_resolve ) {
            labelled_tree_process_edge_fully_resolved_(
                result, *edge_list[ edge_i ], place_map[ edge_i ], name_prefix
            );
        } else {
            labelled_tree_process_edge_multifurcating_(
                result, *edge_list[ edge_i ], place_map[ edge_i ], name_prefix
            );
        }
    }
 
    return result;
}
 
} // namespace placement
} // namespace genesis