A toolkit for working with phylogenetic data.
v0.24.0
tree/function/manipulation.cpp
Go to the documentation of this file.
1 /*
2  Genesis - A toolkit for working with phylogenetic data.
3  Copyright (C) 2014-2019 Lucas Czech 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 
37 #include "genesis/tree/tree.hpp"
39 
42 
43 #include <algorithm>
44 #include <cassert>
45 #include <vector>
46 
47 namespace genesis {
48 namespace tree {
49 
50 // =================================================================================================
51 // Minimal Tree
52 // =================================================================================================
53 
55 {
56  Tree result;
57 
58  // Get container access shortcuts.
59  auto& links = result.expose_link_container();
60  auto& nodes = result.expose_node_container();
61  auto& edges = result.expose_edge_container();
62 
63  // Create all needed elements...
64  auto link1_u = utils::make_unique< TreeLink >();
65  auto link2_u = utils::make_unique< TreeLink >();
66  auto node1_u = utils::make_unique< TreeNode >();
67  auto node2_u = utils::make_unique< TreeNode >();
68  auto edge0_u = utils::make_unique< TreeEdge >();
69 
70  // ... and connect them with each other.
71  link1_u->reset_index( 0 );
72  link1_u->reset_next( link1_u.get() );
73  link1_u->reset_outer( link2_u.get() );
74  link1_u->reset_node( node1_u.get() );
75  link1_u->reset_edge( edge0_u.get() );
76 
77  link2_u->reset_index( 1 );
78  link2_u->reset_next( link2_u.get() );
79  link2_u->reset_outer( link1_u.get() );
80  link2_u->reset_node( node2_u.get() );
81  link2_u->reset_edge( edge0_u.get() );
82 
83  node1_u->reset_index( 0 );
84  node1_u->reset_primary_link( link1_u.get() );
85 
86  node2_u->reset_index( 1 );
87  node2_u->reset_primary_link( link2_u.get() );
88 
89  edge0_u->reset_index( 0 );
90  edge0_u->reset_primary_link( link1_u.get() );
91  edge0_u->reset_secondary_link( link2_u.get() );
92 
93  // Finally, move all elements to the tree and reset the tree root link.
94  links.push_back( std::move( link1_u ));
95  links.push_back( std::move( link2_u ));
96  nodes.push_back( std::move( node1_u ));
97  nodes.push_back( std::move( node2_u ));
98  edges.push_back( std::move( edge0_u ));
99  result.reset_root_link( links.front().get() );
100 
101  return result;
102 }
103 
104 // =================================================================================================
105 // Add Nodes
106 // =================================================================================================
107 
108 TreeNode& add_new_node( Tree& tree, TreeNode& target_node )
109 {
110  // Basic check.
111  if( ! belongs_to( target_node, tree ) ) {
112  throw std::runtime_error(
113  "Cannot add Node to a Tree where the given Node is not part of the Tree."
114  );
115  }
116 
117  // Get container access shortcuts.
118  auto& links = tree.expose_link_container();
119  auto& nodes = tree.expose_node_container();
120  auto& edges = tree.expose_edge_container();
121 
122  // Create all new elements that we need:
123  // 1. A link that gets added to the given node and connects it to the new node.
124  // 2. The link that belongs to the newly created leaf node.
125  // 3. The newly created node itself.
126  // 4. The edge between the two nodes.
127  auto con_link_u = utils::make_unique< TreeLink >();
128  auto end_link_u = utils::make_unique< TreeLink >();
129  auto end_node_u = utils::make_unique< TreeNode >();
130  auto con_edge_u = utils::make_unique< TreeEdge >();
131 
132  // Get the pointers, for ease of use.
133  auto const con_link = con_link_u.get();
134  auto const end_link = end_link_u.get();
135  auto const end_node = end_node_u.get();
136  auto const con_edge = con_edge_u.get();
137 
138  // Populate the link at the given node.
139  con_link->reset_index( links.size() );
140  con_link->reset_node( &target_node );
141  con_link->reset_edge( con_edge );
142  con_link->reset_outer( end_link );
143 
144  // Find the last link of the given node (in traversal order around the node).
145  auto const up_link = &target_node.primary_link();
146  auto last_link = up_link;
147  while( &last_link->next() != up_link ) {
148  last_link = &last_link->next();
149  }
150 
151  // Now insert the new link in between up_link and last_link.
152  assert( &last_link->next() == up_link );
153  last_link->reset_next( con_link );
154  con_link->reset_next( up_link );
155  assert( &last_link->next() == con_link );
156  assert( &con_link->next() == up_link );
157 
158  // Add the link to the tree. This changes the size of the links vector, so the next call to
159  // links.size() will give a new, proper value for the other link.
160  links.push_back( std::move( con_link_u ));
161 
162  // Populate the link at the new node and add it to the tree.
163  end_link->reset_index( links.size() );
164  end_link->reset_node( end_node );
165  end_link->reset_edge( con_edge );
166  end_link->reset_next( end_link );
167  end_link->reset_outer( con_link );
168  links.push_back( std::move( end_link_u ));
169 
170  // Populate the new node and add it to the tree.
171  end_node->reset_index( nodes.size() );
172  end_node->reset_primary_link( end_link );
173  end_node->reset_data( target_node.data_ptr()->recreate() );
174  nodes.push_back( std::move( end_node_u ));
175 
176  // Populate the new edge and add it to the tree.
177  con_edge->reset_index( edges.size() );
178  con_edge->reset_primary_link( con_link );
179  con_edge->reset_secondary_link( end_link );
180  con_edge->reset_data( target_node.primary_link().edge().data_ptr()->recreate() );
181  edges.push_back( std::move( con_edge_u ));
182 
183  // Return the new node. We just moved the node uniq ptr, but not the pointee, so this is valid.
184  return *end_node;
185 }
186 
188  Tree& tree,
189  TreeEdge& target_edge,
190  std::function<void( TreeEdge& target_edge, TreeEdge& new_edge )> adjust_edges
191 ) {
192  // Basic check.
193  if( ! belongs_to( target_edge, tree ) ) {
194  throw std::runtime_error(
195  "Cannot add Node to Tree on an Edge that is not part of the Tree."
196  );
197  }
198 
199  // Get container access shortcuts.
200  auto& links = tree.expose_link_container();
201  auto& nodes = tree.expose_node_container();
202  auto& edges = tree.expose_edge_container();
203 
204  // Create all new elements that we need:
205  // * Two links that build a new node in the middle of the target edge.
206  // * The new node in the middle of the target edge.
207  // * A new edge that connects to the secondary end of the target edge.
208  auto pri_link_u = utils::make_unique< TreeLink >();
209  auto sec_link_u = utils::make_unique< TreeLink >();
210  auto mid_node_u = utils::make_unique< TreeNode >();
211  auto sec_edge_u = utils::make_unique< TreeEdge >();
212 
213  // Get the pointers, for ease of use.
214  auto const pri_link = pri_link_u.get();
215  auto const sec_link = sec_link_u.get();
216  auto const mid_node = mid_node_u.get();
217  auto const sec_edge = sec_edge_u.get();
218 
219  // Populate the link towards the primary end of the target edge and add it to the tree.
220  pri_link->reset_index( links.size() );
221  pri_link->reset_next( sec_link );
222  pri_link->reset_outer( &target_edge.primary_link() );
223  pri_link->reset_node( mid_node );
224  pri_link->reset_edge( &target_edge );
225  links.push_back( std::move( pri_link_u ));
226 
227  // Populate the link towards the secondary end of the target edge and add it to the tree.
228  sec_link->reset_index( links.size() );
229  sec_link->reset_next( pri_link );
230  sec_link->reset_outer( &target_edge.secondary_link() );
231  sec_link->reset_node( mid_node );
232  sec_link->reset_edge( sec_edge );
233  links.push_back( std::move( sec_link_u ));
234 
235  // Populate the new node in the middle of the target edge and add it to the tree.
236  mid_node->reset_index( nodes.size() );
237  mid_node->reset_primary_link( pri_link );
238  mid_node->reset_data( target_edge.primary_node().data_ptr()->recreate() );
239  nodes.push_back( std::move( mid_node_u ));
240 
241  // Populate the edge at the secondary end of the target edge and add it to the tree.
242  sec_edge->reset_index( edges.size() );
243  sec_edge->reset_primary_link( sec_link );
244  sec_edge->reset_secondary_link( &target_edge.secondary_link() );
245  sec_edge->reset_data( target_edge.data_ptr()->recreate() );
246  edges.push_back( std::move( sec_edge_u ));
247 
248  // Finally adjust the existing tree elements.
249  target_edge.primary_link().reset_outer( pri_link );
250  target_edge.secondary_link().reset_outer( sec_link );
251  target_edge.secondary_link().reset_edge( sec_edge );
252  target_edge.reset_secondary_link( pri_link );
253 
254  // If we have a transform function, call it.
255  if( adjust_edges ) {
256  adjust_edges( target_edge, *sec_edge );
257  }
258 
259  return *mid_node;
260 }
261 
263  Tree& tree,
264  TreeEdge& target_edge,
265  std::function<void( TreeEdge& target_edge, TreeEdge& new_edge )> adjust_edges
266 ) {
267  // First add a node that splits the edge, and then a new leaf node to this one.
268  auto& mid_node = add_new_node( tree, target_edge, adjust_edges );
269  return add_new_node( tree, mid_node );
270 }
271 
272 // =================================================================================================
273 // Delete Nodes
274 // =================================================================================================
275 
276 void delete_node( Tree& tree, TreeNode& target_node )
277 {
278  // Basic check.
279  if( ! belongs_to( target_node, tree ) ) {
280  throw std::runtime_error(
281  "Cannot delete Node from a Tree that is not part of the Tree."
282  );
283  }
284 
285  auto const deg = degree( target_node );
286  if( deg == 1 ) {
287  delete_leaf_node( tree, target_node );
288  } else if( deg == 2 ) {
289  delete_linear_node( tree, target_node );
290  } else {
291  delete_subtree( tree, Subtree{ target_node } );
292  }
293 }
294 
295 void delete_leaf_node( Tree& tree, TreeNode& target_node )
296 {
297  // Basic check.
298  if( ! belongs_to( target_node, tree ) ) {
299  throw std::runtime_error(
300  "Cannot delete node from a tree that is not part of the tree."
301  );
302  }
303  if( degree( target_node ) != 1 ) {
304  throw std::runtime_error(
305  "Cannot delete leaf node if the node is not actually a leaf."
306  );
307  }
308  if( tree.link_count() <= 2 ) {
309  assert( tree.node_count() == 2 && tree.edge_count() == 1 );
310  throw std::runtime_error(
311  "Cannot delete leaf node from a minmal tree that only contains two nodes."
312  );
313  }
314 
315  // Node has to be a leaf.
316  assert( &target_node.link().next() == &target_node.link() );
317 
318  // Get container access shortcuts.
319  auto& links = tree.expose_link_container();
320  auto& nodes = tree.expose_node_container();
321  auto& edges = tree.expose_edge_container();
322 
323  // If the node to be deleted is the root, we need to reset to its adjacent node.
324  // We need to check this now, and reset later, because the link indices will change...
325  // Also, if the root is the node where the subtree is attached, we might need to reset
326  // the root link, because the attachement link might get deleted.
327  // For this, we use the first link of the node where the node is attached.
328  auto root_link = &tree.root_link();
329  if( tree.root_node().index() == target_node.index() || &target_node.link().outer() == &tree.root_link() ) {
330  root_link = &target_node.link().outer().next();
331 
332  // Also, if we reset the root link, we need to update the primary link of the new root node.
333  // Do some checks for this. (The commented line is an alternative to resetting it.
334  // We now do this later, at the end of the function.)
335  assert( &target_node.link().outer().node().link() == &target_node.link().outer() );
336  assert( &root_link->node() == &target_node.link().outer().node() );
337  // target_node.link().outer().node().reset_primary_link( root_link );
338  }
339 
340  // Delete the edge and adjust indices of other edges as needed.
341  // We do this first, so that the link to the edge is still valid.
342  auto const edge_idx = target_node.link().edge().index();
343  edges.erase( edges.begin() + edge_idx );
344  for( size_t i = edge_idx; i < tree.edge_count(); ++i ) {
345  assert( tree.edge_at(i).index() == i + 1 );
346  tree.edge_at(i).reset_index(i);
347  }
348 
349  // Make the node that the target is attached to forget about this node by skipping the link.
350  auto& attach_link = target_node.link().outer();
351  assert( &attach_link.edge() == &target_node.link().edge() );
352  assert( &attach_link.outer() == &target_node.link() );
353  auto link_ptr = &( attach_link.next() );
354  while( &link_ptr->next() != &attach_link ) {
355  link_ptr = &link_ptr->next();
356  }
357  assert( &link_ptr->next() == &attach_link );
358  assert( &link_ptr->next().next() == &attach_link.next() );
359  link_ptr->reset_next( &link_ptr->next().next() );
360 
361  // Delete both links, and adjust the indices of the other links as needed.
362  // Deleting two elemetns is tricky, as the indices shift...
363  auto const minmax_link_idx = utils::minmax_value( attach_link.index(), attach_link.outer().index() );
364  links.erase( links.begin() + minmax_link_idx.first );
365  links.erase( links.begin() + minmax_link_idx.second - 1 );
366  for( size_t i = minmax_link_idx.first; i < links.size(); ++i ) {
367  assert( tree.link_at(i).index() == i + 1 || tree.link_at(i).index() == i + 2 );
368  tree.link_at(i).reset_index(i);
369  }
370 
371  // Finally, delete the node and adjust indices of other nodes as needed.
372  auto const node_idx = target_node.index();
373  nodes.erase( nodes.begin() + node_idx );
374  for( size_t i = node_idx; i < tree.node_count(); ++i ) {
375  assert( tree.node_at(i).index() == i + 1 );
376  tree.node_at(i).reset_index(i);
377  }
378 
379  // Lastly, reset the root link (or set it to what it was before...)
380  tree.reset_root_link( root_link );
381  tree.root_link().node().reset_primary_link( root_link );
382 }
383 
385  Tree& tree,
386  TreeNode& target_node,
387  std::function<void( TreeEdge& remaining_edge, TreeEdge& deleted_edge )> adjust_edges
388 ) {
389  // Basic check.
390  if( ! belongs_to( target_node, tree ) ) {
391  throw std::runtime_error(
392  "Cannot delete Node from a Tree that is not part of the Tree."
393  );
394  }
395  if( degree( target_node ) != 2 ) {
396  throw std::runtime_error(
397  "Cannot delete linear Node if the Node is not actually linear (degree 2)."
398  );
399  }
400 
401  // Before we do anything, call the adjust function.
402  if( adjust_edges ) {
403  adjust_edges( target_node.primary_link().edge(), target_node.primary_link().next().edge() );
404  }
405 
406  // If the node to be deleted is the root, we need to reset to its adjacent node.
407  // We need to check this now, and reset later, because the node indices will change...
408  auto root_link = &tree.root_link();
409  if( tree.root_node().index() == target_node.index() ) {
410  root_link = &target_node.link().outer();
411  }
412 
413  // Get container access shortcuts.
414  auto& links = tree.expose_link_container();
415  auto& nodes = tree.expose_node_container();
416  auto& edges = tree.expose_edge_container();
417 
418  // Adjust the links of the two adjacent nodes to point to each other.
419  auto& pr_link = target_node.link();
420  auto& adj_link_p = pr_link.outer();
421  auto& adj_link_d = pr_link.next().outer();
422  assert( &adj_link_p.outer().node() == &adj_link_d.outer().node() );
423  adj_link_p.reset_outer( &adj_link_d );
424  adj_link_d.reset_outer( &adj_link_p );
425 
426  // Adjust the pointers to and from the remaining edge.
427  // We are resetting the primary adjacent edge here. This also takes care of setting it to
428  // another node in case that the deleted node is the root.
429  // The first assertion needs a special case for this (if the deletion node is the root).
430  assert( &adj_link_p.edge().secondary_node() == &target_node || root_link != &tree.root_link() );
431  assert( &adj_link_p.edge() == &target_node.link().edge() );
432  assert( &adj_link_d.edge() == &target_node.link().next().edge() );
433  assert( &adj_link_d.edge().primary_node() == &target_node );
434  adj_link_p.edge().reset_primary_link( &adj_link_p );
435  adj_link_p.edge().reset_secondary_link( &adj_link_d );
436  adj_link_d.reset_edge( &target_node.primary_link().edge() );
437 
438  // Delete the edge and adjust indices of other edges as needed.
439  auto const edge_idx = pr_link.next().edge().index();
440  edges.erase( edges.begin() + edge_idx );
441  for( size_t i = edge_idx; i < tree.edge_count(); ++i ) {
442  assert( tree.edge_at(i).index() == i + 1 );
443  tree.edge_at(i).reset_index(i);
444  }
445 
446  // Delete both links of the node, and adjust the indices of the other links as needed.
447  // Deleting two elemetns is tricky, as the indices shift...
448  auto const minmax_link_idx = utils::minmax_value( pr_link.index(), pr_link.next().index() );
449  links.erase( links.begin() + minmax_link_idx.first );
450  links.erase( links.begin() + minmax_link_idx.second - 1 );
451  for( size_t i = minmax_link_idx.first; i < links.size(); ++i ) {
452  assert( tree.link_at(i).index() == i + 1 || tree.link_at(i).index() == i + 2 );
453  tree.link_at(i).reset_index(i);
454  }
455 
456  // Finally, delete the node and adjust indices of other nodes as needed.
457  auto const node_idx = target_node.index();
458  nodes.erase( nodes.begin() + node_idx );
459  for( size_t i = node_idx; i < tree.node_count(); ++i ) {
460  assert( tree.node_at(i).index() == i + 1 );
461  tree.node_at(i).reset_index(i);
462  }
463 
464  // Lastly, reset the root link (or set it to what it was before...)
465  tree.reset_root_link( root_link );
466 }
467 
474 template<class T>
475 static void delete_from_tree_container_( T& old_elems, std::vector<size_t>& del_elems )
476 {
477  // Sort the deletion list, so that we do not need to seek in it. There should be no duplicats.
478  std::sort( del_elems.begin(), del_elems.end() );
479  assert( std::adjacent_find( del_elems.begin(), del_elems.end()) == del_elems.end() );
480 
481  // We never delete the whole tree!
482  assert( del_elems.size() < old_elems.size() );
483 
484  // Create an empty container and move the needed elements there.
485  // This is just so much easier than trying to delete within a container
486  // while iterating it, adjusting the internal indices of the elemetns, etc.
487  auto new_elems = T{};
488  new_elems.reserve( old_elems.size() - del_elems.size() );
489 
490  // Move elements and adjust indices.
491  size_t del_elems_idx = 0;
492  for( size_t i = 0; i < old_elems.size(); ++i ) {
493  assert( i == old_elems[i]->index() );
494  if( del_elems_idx < del_elems.size() && del_elems[del_elems_idx] == i ) {
495 
496  // If the current index is in the deletion list, don't move the element,
497  // and go to the next element of the deletion list.
498  ++del_elems_idx;
499  } else {
500 
501  // If the current index is not in the deletion list, move it to the new container.
502  assert( old_elems[i]->index() == i );
503  assert( new_elems.size() == i - del_elems_idx );
504  old_elems[i]->reset_index( i - del_elems_idx );
505  new_elems.emplace_back( std::move( old_elems[i] ));
506  assert( new_elems.back()->index() == new_elems.size() - 1 );
507  }
508  }
509 
510  // Swap, which will put the elements that we want to delete in the new container,
511  // which will destroy them when going out of scope.
512  assert( new_elems.size() + del_elems.size() == old_elems.size() );
513  std::swap( new_elems, old_elems );
514  // old_elems = std::move( new_elems );
515 }
516 
517 void delete_subtree( Tree& tree, Subtree const& subtree )
518 {
519  // Basic check.
520  if( ! belongs_to( subtree, tree ) ) {
521  throw std::runtime_error(
522  "Cannot delete Subtree from a Tree that is not part of the Tree."
523  );
524  }
525  if( is_leaf( subtree.link().outer() )) {
526  throw std::runtime_error(
527  "Cannot delete Subtree from Tree if this leads to a singluar Node."
528  );
529  }
530 
531  // Make lists of all indices that need to be deleted.
532  // This already includes the edge towards the rest of the tree, as well as the link
533  // at the attachment node!
534  // Also, find out if whe are deleting the root. If so, we need to reset it later.
535  std::vector<size_t> del_nodes;
536  std::vector<size_t> del_edges;
537  std::vector<size_t> del_links;
538  bool contains_root = false;
539  for( auto it : preorder( subtree )) {
540  del_nodes.push_back( it.node().index() );
541  del_edges.push_back( it.edge().index() );
542  del_links.push_back( it.link().index() );
543  del_links.push_back( it.link().outer().index() );
544 
545  if( is_root( it.node() )) {
546  contains_root = true;
547  }
548  }
549 
550  // If we are about to delete the root, store the new link for later.
551  // Also, if the root is the node where the subtree is attached, we might need to reset
552  // the root link, because the attachement link might get deleted.
553  // For this, we use the first link of the node where the subtree is attached.
554  // This cannot be a leaf node (as we already checked this earlier), otherwise we'd be setting
555  // it to the attachement link itself, because if this was the case, it's next() is itself.
556  // If not, we just store the current one, so that setting it later doesn't change anything.
557  assert( ! is_leaf( subtree.link().outer() ));
558  auto root_link = &tree.root_link();
559  if( contains_root || &subtree.link().outer() == &tree.root_link() ) {
560  root_link = &tree.link_at( subtree.link().outer().next().index() );
561  }
562 
563  // Get the link that points to the attachment link
564  // (the one on the remaining node that will be deleted).
565  auto& attach_link = subtree.link().outer();
566  assert( &attach_link.edge() == &subtree.link().edge() );
567  assert( &attach_link.outer() == &subtree.link() );
568  auto link_ptr = &( attach_link.next() );
569  while( &link_ptr->next() != &attach_link ) {
570  link_ptr = &link_ptr->next();
571  }
572  assert( &link_ptr->next() == &attach_link );
573  assert( &link_ptr->next().next() == &attach_link.next() );
574  auto& link_nc = tree.link_at( link_ptr->index() );
575 
576  // If the primary link of the attachment node is about to be deleted, reset it to another link.
577  if( &link_nc.node().primary_link() == &link_nc.next() ) {
578  assert( &attach_link.node().primary_link() == &attach_link );
579  link_nc.node().reset_primary_link( &link_nc.next().next() );
580  }
581 
582  // Make the node that the target is attached to forget about this node by skipping the link.
583  // Because the subtree always stores a const link, we need to get a non-const from the tree.
584  link_nc.reset_next( &link_nc.next().next() );
585 
586  // Delete the elements at the given indices.
590 
591  // Reset the root if neeeded. This uses the new index of the pointee.
592  tree.reset_root_link( root_link );
593  tree.root_link().node().reset_primary_link( root_link );
594 }
595 
596 // void delete_edge(
597 // Tree& tree,
598 // TreeEdge& target_edge,
599 // std::function<void( TreeEdge& remaining_node, TreeEdge& deleted_node )> adjust_nodes
600 // ) {
601 //
602 // }
603 
604 // =================================================================================================
605 // Rerooting
606 // =================================================================================================
607 
609  Tree& tree,
610  TreeEdge& target_edge,
611  std::function<void( TreeEdge& target_edge, TreeEdge& new_edge )> adjust_edges
612 ) {
613  if( is_rooted( tree )) {
614  throw std::runtime_error( "Cannot make a Tree rooted if it is already rooted." );
615  }
616 
617  auto& mid_node = add_new_node( tree, target_edge, adjust_edges );
618  change_rooting( tree, mid_node );
619  return mid_node;
620 }
621 
623  Tree& tree,
624  std::function<void( TreeEdge& target_edge, TreeEdge& new_edge )> adjust_edges
625 ) {
626  return make_rooted( tree, tree.root_link().edge(), adjust_edges );
627 }
628 
630  Tree& tree,
631  std::function<void( TreeEdge& remaining_edge, TreeEdge& deleted_edge )> adjust_edges
632 ) {
633  if( ! is_rooted( tree )) {
634  throw std::runtime_error( "Cannot make a Tree unrooted if it is already unrooted." );
635  }
636 
637  delete_linear_node( tree, tree.root_node(), adjust_edges );
638 }
639 
640 void change_rooting( Tree& tree, TreeLink& at_link )
641 {
642  if( ! belongs_to( at_link, tree ) ) {
643  throw std::runtime_error( "Cannot reroot Tree on a Link that is not part of the Tree." );
644  }
645 
646  // We store the old root node, becasue we will change internals of the tree, so that
647  // node().is_root() won't work while this function is running.
648  TreeNode* old_root = &tree.root_node();
649 
650  // Pointer to the primary link of the target node.
651  TreeLink* cur_link = &tree.link_at( at_link.index() ).node().primary_link();
652 
653  // Set new root index and node link of the new root.
654  assert( &at_link == &tree.link_at( at_link.index() ));
655  tree.reset_root_link( &at_link );
656  at_link.node().reset_primary_link( &tree.link_at( at_link.index() ));
657 
658  // Walk the path from the new root to the old, and change all pointers of the edges and nodes
659  // on that path so that they point towards the new root.
660  while( &cur_link->node() != old_root ) {
661 
662  // Assert that the primary direction is correct: Is the current link is at the secondary
663  // end of its edge?
664  assert( cur_link == &cur_link->edge().secondary_link() );
665 
666  // Swap the edge's links, so that they point towards the new root.
667  auto link_p_tmp = &cur_link->edge().primary_link();
668  auto link_s_tmp = &cur_link->edge().secondary_link();
669  cur_link->edge().reset_primary_link( link_s_tmp );
670  cur_link->edge().reset_secondary_link( link_p_tmp );
671  // std::swap( cur_link->edge_->link_p_, cur_link->edge_->link_s_ );
672 
673  // Assert that this worked.
674  assert( cur_link == &cur_link->edge().primary_link() );
675  assert( &cur_link->outer() == &cur_link->edge().secondary_link() );
676 
677  // Store the link of the next node that points towards the root.
678  // We need it, because we will change this upwards link of the next node now.
679  auto to_root_link = &cur_link->outer().node().primary_link();
680 
681  // Change the main link of the next node so that it points towards the new root.
682  cur_link->outer().node().reset_primary_link( &cur_link->outer() );
683 
684  // Move one node towards the root.
685  cur_link = to_root_link;
686  }
687 }
688 
689 void change_rooting( Tree& tree, TreeNode& at_node )
690 {
691  if( ! belongs_to( at_node, tree ) ) {
692  throw std::runtime_error( "Cannot reroot Tree on a Node that is not part of the Tree." );
693  }
694  change_rooting( tree, at_node.link() );
695 }
696 
697 // =================================================================================================
698 // Ladderize
699 // =================================================================================================
700 
701 void ladderize( Tree& tree, LadderizeOrder order )
702 {
703  // For each node, get how many nodes its subtree (away from the root) has.
704  // We use this quantity to sort each node's links.
705  auto sub_sizes = subtree_sizes( tree );
706 
707  // Ladderize all nodes
708  for( auto& node : tree.nodes() ) {
709 
710  // No need to ladderize a leaf. It would still work, but we can use this as a speedup.
711  if( is_leaf( node )) {
712  continue;
713  }
714 
715  // Get the sizes of the children/subtrees of this node.
716  std::vector<size_t> child_sizes;
717  std::vector<TreeLink*> child_links;
718  for( auto link_it : node_links( node ) ) {
719 
720  // Don't treat the link towards the root; we only want to sort the subtree.
721  // Assert that the first iteration is actually this link towards the root.
722  if( link_it.is_first_iteration() ) {
723  assert( &link_it.link() == &node.primary_link() );
724  continue;
725  }
726 
727  child_sizes.push_back( sub_sizes[ link_it.link().outer().node().index() ] );
728  child_links.push_back( &link_it.link() );
729  }
730 
731  // Sort the sizes. We use stable sort in order to not change the order of equal sized subtrees.
732  auto child_order = ( order == LadderizeOrder::kSmallFirst
733  ? utils::stable_sort_indices( child_sizes.begin(), child_sizes.end(), std::less<size_t>() )
734  : utils::stable_sort_indices( child_sizes.begin(), child_sizes.end(), std::greater<size_t>() )
735  );
736 
737  // The number of indices needs to be the rank of the node (number of immediate children).
738  assert( child_order.size() == child_sizes.size() );
739  assert( child_order.size() == child_links.size() );
740  assert( child_order.size() == degree( node ) - 1 );
741 
742  // Change all next links of the node so that they reflect the subtree size order.
743  auto cur_link = &node.primary_link();
744  for( auto child_order_i : child_order ) {
745 
746  // We use this assertion to ensure that each link is only processed once.
747  // At the end of this loop, we set it to nullptr, so a second encounter would fail.
748  assert( child_links[ child_order_i ] );
749 
750  // Set the link's next link and move on.
751  cur_link->reset_next( child_links[ child_order_i ] );
752  cur_link = child_links[ child_order_i ];
753 
754  // Set the link in the vector to null, so that the above assert can check that we
755  // never process it again.
756  child_links[ child_order_i ] = nullptr;
757  }
758 
759  // We now need to set the next pointer of the last link of the node so that it points
760  // back to the original starting node (the one towards the root).
761  cur_link->reset_next( &node.primary_link() );
762 
763  // Finally, assert that we processed all links. If so, all of them are null by now.
764  for( auto const& cl : child_links ) {
765  (void) cl;
766  assert( !cl );
767  }
768  }
769 }
770 
771 } // namespace tree
772 } // namespace genesis
void delete_node(Tree &tree, TreeNode &target_node)
Delete a TreeNode from a Tree.
utils::Range< IteratorNodeLinks< true > > node_links(ElementType const &element)
Definition: node_links.hpp:186
TreeLink & secondary_link()
Return the TreeLink of this TreeEdge that points away from the root.
Definition: edge.hpp:130
bool is_rooted(Tree const &tree)
Return whether the Tree is rooted, that is, whether the root node has two neighbors.
Provides some valuable algorithms that are not part of the C++ 11 STL.
TreeNode & reset_primary_link(TreeLink *val)
Reset the internal pointer to the TreeLink of this TreeNode.
Definition: node.hpp:276
void ladderize(Tree &tree, LadderizeOrder order)
Ladderize a Tree, that is, order its subtrees by size.
virtual std::unique_ptr< BaseEdgeData > recreate() const
Polymorphically create a default-constructed instance of this class with the same derived type as it ...
Definition: edge_data.hpp:126
bool is_root(TreeLink const &link)
Return whether the link belongs to the root node of its Tree.
void make_unrooted(Tree &tree, std::function< void(TreeEdge &remaining_edge, TreeEdge &deleted_edge)> adjust_edges)
Unroot a Tree.
void swap(SequenceSet &lhs, SequenceSet &rhs)
void delete_subtree(Tree &tree, Subtree const &subtree)
Delete a complete Subtree from a Tree.
Tree operator functions.
TreeLink & primary_link()
Return the TreeLink that points towards the root.
Definition: node.hpp:110
Tree minimal_tree_topology()
Create a minimal Tree that can be used with manipulation functions such as add_new_node() or add_new_...
size_t edge_count() const
Return the number of TreeEdges of the Tree.
Definition: tree/tree.hpp:272
size_t index() const
Return the index of this Edge.
Definition: edge.hpp:106
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
Reference to a subtree of a Tree.
Definition: subtree.hpp:69
Provides some valuable additions to STD.
TreeNode & make_rooted(Tree &tree, TreeEdge &target_edge, std::function< void(TreeEdge &target_edge, TreeEdge &new_edge)> adjust_edges)
Root a Tree at a given TreeEdge.
TreeNode & primary_node()
Return the TreeNode of this TreeEdge that points towards the root.
Definition: edge.hpp:146
static void delete_from_tree_container_(T &old_elems, std::vector< size_t > &del_elems)
Local helper function template that takes one of the tree element containers and deletes all elements...
BaseEdgeData * data_ptr()
Return a pointer to the data.
Definition: edge.hpp:246
utils::Range< IteratorNodes > nodes()
Definition: tree/tree.hpp:418
Class for representing phylogenetic trees.
Definition: tree/tree.hpp:97
size_t link_count() const
Return the number of TreeLinks of the Tree.
Definition: tree/tree.hpp:256
TreeEdge & reset_index(size_t val)
Reset the internal index of this TreeEdge.
Definition: edge.hpp:276
TreeLink & link()
Return the TreeLink that points towards the root.
Definition: node.hpp:129
Tree & reset_root_link(TreeLink *root_link)
Reset the link that is considered to be the root of the Tree.
Definition: tree/tree.cpp:184
virtual std::unique_ptr< BaseNodeData > recreate() const
Polymorphically create a default-constructed instance of this class with the same derived type as it ...
Definition: node_data.hpp:126
TreeNode & node_at(size_t index)
Return the TreeNode at a certain index.
Definition: tree/tree.hpp:220
TreeLink & root_link()
Return the TreeLink at the current root of the Tree.
Definition: tree/tree.hpp:284
LinkContainerType & expose_link_container()
Get the container that stores all TreeLinks of the Tree.
Definition: tree/tree.cpp:192
TreeEdge & reset_primary_link(TreeLink *val)
Reset the internal pointer to the primary TreeLink of this TreeEdge.
Definition: edge.hpp:290
TreeNode & add_new_leaf_node(Tree &tree, TreeEdge &target_edge, std::function< void(TreeEdge &target_edge, TreeEdge &new_edge)> adjust_edges)
Add a new Node as a leaf to an existing Edge, by also adding a new Node in the middle of that Edge...
TreeEdge & edge_at(size_t index)
Return the TreeEdge at a certain index.
Definition: tree/tree.hpp:238
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...
TreeLink const & link() const
Get the TreeLink that separates the subtree from the rest of the tree.
Definition: subtree.hpp:125
TreeNode & add_new_node(Tree &tree, TreeNode &target_node)
Add a new Node as a leaf to an existing Node.
std::vector< size_t > stable_sort_indices(RandomAccessIterator first, RandomAccessIterator last, Comparator comparator)
Get the indices to the stable sorted order of the given range.
Definition: algorithm.hpp:239
size_t node_count() const
Return the number of TreeNodes of the Tree.
Definition: tree/tree.hpp:264
utils::Range< IteratorPreorder< true > > preorder(ElementType const &element)
TreeLink & link_at(size_t index)
Return the TreeLink at a certain index.
Definition: tree/tree.hpp:202
TreeNode & root_node()
Return the TreeNode at the current root of the Tree.
Definition: tree/tree.hpp:300
void delete_leaf_node(Tree &tree, TreeNode &target_node)
Delete a leaf TreeNode.
std::pair< T, T > minmax_value(T const &a, T const &b)
Returns the lowest and the greatest of the given values, by value.
Definition: algorithm.hpp:128
Header of Tree class.
size_t index() const
Return the index of this Node.
Definition: node.hpp:102
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...
void change_rooting(Tree &tree, TreeLink &at_link)
"Reroot" the Tree at the given TreeLink.
void delete_linear_node(Tree &tree, TreeNode &target_node, std::function< void(TreeEdge &remaining_edge, TreeEdge &deleted_edge)> adjust_edges)
Delete a "linear" TreeNode from a Tree, that is, a node with two neighbours.
TreeNode & reset_index(size_t val)
Reset the internal index of this TreeNode.
Definition: node.hpp:262
EdgeContainerType & expose_edge_container()
Get the container that stores all TreeEdges of the Tree.
Definition: tree/tree.cpp:202
BaseNodeData * data_ptr()
Return a pointer to the data.
Definition: node.hpp:232
bool is_leaf(TreeLink const &link)
Return true iff the node of the given link is a leaf node.
TreeEdge & reset_secondary_link(TreeLink *val)
Reset the internal pointer to the secondary TreeLink of this TreeEdge.
Definition: edge.hpp:304
NodeContainerType & expose_node_container()
Get the container that stores all TreeNodes of the Tree.
Definition: tree/tree.cpp:197
std::vector< size_t > subtree_sizes(Tree const &tree, TreeNode const &node)
Calculate the sizes of all subtrees as seen from the given TreeNode.