A toolkit for working with phylogenetic data.
v0.24.0
placement/function/manipulation.cpp
Go to the documentation of this file.
1 /*
2  Genesis - A toolkit for working with phylogenetic data.
3  Copyright (C) 2014-2018 Lucas Czech, Pierre Barbera and HITS gGmbH
4 
5  This program is free software: you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation, either version 3 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 
18  Contact:
19  Lucas Czech <lucas.czech@h-its.org>
20  Exelixis Lab, Heidelberg Institute for Theoretical Studies
21  Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
22 */
23 
32 
36 
41 
42 #include <algorithm>
43 #include <cassert>
44 #include <stdexcept>
45 #include <unordered_set>
46 #include <vector>
47 
48 namespace genesis {
49 namespace placement {
50 
51 void make_rooted( Sample& sample, PlacementTreeEdge& target_edge )
52 {
53  auto& tree = sample.tree();
54 
55  if( tree::is_rooted( tree )) {
56  throw std::invalid_argument{ "Cannot root a PlacementTree that is already rooted." };
57  }
58  if( not tree::belongs_to( tree, target_edge )) {
59  throw std::invalid_argument{ "The given edge does not belong to the tree." };
60  }
61 
62  // Get old root
63  auto& old_root = tree.root_node();
64 
65  // Make rooted: add root to edge. Make sure that the new node points
66  auto& new_node = make_rooted( tree, target_edge );
67  auto& new_edge = target_edge.primary_link().next().edge();
68  assert( &new_node.primary_link() == &target_edge.primary_link() );
69  assert( degree( new_node ) == 2 );
70  (void) new_node;
71 
72  // This creates a new TreeNode, with the original `edge` being one of the two edges adjacent to the new node.
73  // the original edge will be the one that is more _towards_ the original root of the tree,
74  // see add_new_node( Tree&, TreeEdge&, ...)
75  // Thus, as we generally want this new root to be more toward the interior of the tree,
76  // we scale the branch length accordingly.
77 
78  // Rescale adjacent branch lengths of newly created root node to 0% and 100%
79  auto const original_branch_length = target_edge.data<PlacementEdgeData>().branch_length;
80  target_edge.data<PlacementEdgeData>().branch_length = 0.0;
81  assert( new_edge.has_data() && new_edge.data<PlacementEdgeData>().branch_length == 0.0 );
82  new_edge.data<PlacementEdgeData>().branch_length = original_branch_length;
83 
84  // Next, we need to identify the edges that had their direction to root changed,
85  // as this is information used in the placements (`distal_length` or `proximal_length`).
86  // Once we know which they are, we can iterate over all placements and adjust those numbers
87  // for all placements associated with those edges.
88 
89  // Iterate over the path between old and new root.
90  // path_to_root() also returns the root_link of the tree, but we can already get the
91  // relevant edge from the link before. So out with it
92  auto path = tree::path_to_root( old_root );
93  assert( ! path.empty() && is_root( *path.back() ));
94  path.pop_back();
95 
96  // Get the set of edges towards the old root.
97  std::unordered_set<PlacementTreeEdge*> edges_to_adjust;
98  for( auto link_ptr : path ) {
99  edges_to_adjust.insert( const_cast<PlacementTreeEdge*>( &link_ptr->edge() ) );
100  }
101 
102  // Look for relevant placements, adjust the distal length
103  for( auto& pq : sample ) {
104  for( auto& p : pq.placements() ) {
105  auto place_edge_ptr = &( p.edge() );
106 
107  // If the placement points to the edge on which we rooted, change to new edge
108  if ( place_edge_ptr == &target_edge ) {
109  p.reset_edge( new_edge );
110 
111  // The target edge is on the path. Check this, and change the current pointer
112  // to the new edge, which is not on the path. Otherwise we'd wrongly
113  // flip the proximal length later...
114  assert( edges_to_adjust.count( place_edge_ptr ) > 0 );
115  place_edge_ptr = &( p.edge() );
116  assert( edges_to_adjust.count( place_edge_ptr ) == 0 );
117  }
118 
119  // The current edge can never be the target edge, because we excluded this case above.
120  assert( place_edge_ptr != &target_edge );
121 
122  // If this placement belongs to one of the relevant edges, adjust the proximal_length.
123  if( edges_to_adjust.count( place_edge_ptr ) > 0 ) {
124  auto full_length = place_edge_ptr->data<PlacementEdgeData>().branch_length;
125  p.proximal_length = full_length - p.proximal_length;
126  }
127  }
128  }
129 
130  // Recalculate the edge nums. As the placements use pointer to their edges to get the edge nums,
131  // no change is needed for the placements themselves.
132  reset_edge_nums( tree );
133 }
134 
135 } // namespace placement
136 } // namespace genesis
bool is_rooted(Tree const &tree)
Return whether the Tree is rooted, that is, whether the root node has two neighbors.
Data class for PlacementTreeEdges. Stores the branch length of the edge, and the edge_num, as defined in the jplace standard.
PlacementTree & tree()
Get the PlacementTree of this Sample.
Definition: sample.cpp:119
bool is_root(TreeLink const &link)
Return whether the link belongs to the root node of its Tree.
Tree operator functions.
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
std::vector< TreeLink const *> path_to_root(TreeNode const &node)
Helper function that finds all TreeLinks between a given TreeNode and the root of the Tree...
utils::Range< IteratorPath< true > > path(ElementType const &start, ElementType const &finish)
Definition: path.hpp:337
void reset_edge_nums(PlacementTree &tree)
Reset all edge nums of a PlacementTree.
void make_rooted(Sample &sample, PlacementTreeEdge &target_edge)
Root the underlying PlacementTree of a Sample at a specified TreeEdge.
bool belongs_to(Tree const &tree, TreeNode const &node)
Return whether the TreeNode belongs to the Tree, i.e., whether it is owned by the Tree...
Manage a set of Pqueries along with the PlacementTree where the PqueryPlacements are placed on...
Definition: sample.hpp:68
TreeNode & root_node()
Return the TreeNode at the current root of the Tree.
Definition: tree/tree.hpp:300
TreeLink & primary_link()
Return the TreeLink of this TreeEdge that points towards the root.
Definition: edge.hpp:114
size_t degree(TreeLink const &link)
Return the degree of the node for a given TreeLink, i.e. how many neighbouring nodes it has...
EdgeDataType & data()
Definition: edge.hpp:217
Header of PqueryPlacement class.