57 bool const fully_resolve,
58 std::string
const& name_prefix
73 std::string
const& name_prefix
86 throw std::runtime_error(
87 "Tree provided for labelled_tree() is not compatible "
88 "with the Tree of the provided Sample."
93 std::vector< tree::TreeEdge* > edge_list;
94 edge_list.reserve( result.edge_count() );
95 for(
auto& edge : result.edges() ) {
96 edge_list.push_back( &edge );
115 auto placement_pairs_per_edge = [](
Sample const& sample ){
116 std::vector< std::vector< PlacementPair >> place_map;
118 for(
auto const& pqry : sample.
pqueries() ) {
121 auto max_it = std::max_element(
122 pqry.placements().begin(),
123 pqry.placements().end(),
125 return lhs.like_weight_ratio < rhs.like_weight_ratio;
131 if( max_it != pqry.placements().end() ) {
132 place_map[ max_it->edge().index() ].emplace_back( &pqry, &*max_it );
150 auto place_map = placement_pairs_per_edge( sample );
151 assert( place_map.size() == edge_list.size() );
157 auto sort_placements = fully_resolve
158 ? []( PlacementPair
const& lhs, PlacementPair
const& rhs ){
159 return lhs.placement->proximal_length < rhs.placement->proximal_length;
161 : []( PlacementPair
const& lhs, PlacementPair
const& rhs ){
162 return lhs.placement->pendant_length < rhs.placement->pendant_length;
165 for(
auto& edge_l : place_map ) {
166 std::sort( edge_l.begin(), edge_l.end(), sort_placements );
174 auto add_lonely_placement = [ &name_prefix ](
177 PlacementPair
const& placement_pair
181 auto& pendant_edge = pendant_node.link().edge();
182 auto& proximal_edge = pendant_edge.primary_link().next().edge();
183 auto& distal_edge = pendant_edge.primary_link().next().next().edge();
186 assert(
degree( pendant_edge.primary_node() ) == 3 );
187 assert(
is_leaf( pendant_edge.secondary_node() ));
195 assert( pend_data.branch_length == 0.0 );
199 double const distal_len = prox_data.
branch_length - placement_pair.placement->proximal_length;
200 pend_data.branch_length = placement_pair.placement->pendant_length;
201 prox_data.
branch_length = placement_pair.placement->proximal_length;
206 assert( placement_pair.pquery->name_size() <= 1 );
207 if( placement_pair.pquery->name_size() == 1 ) {
209 pendant_node_data.
name = name_prefix + placement_pair.pquery->name_at(0).name;
218 auto process_edge_fully_resolved = [ &name_prefix ] (
221 std::vector<PlacementPair>
const& placement_pairs
227 double used_length = 0.0;
234 auto insertion_edge = &edge;
237 for(
auto const& placement_pair : placement_pairs ) {
239 for(
auto const& pquery_name : placement_pair.pquery->names() ) {
243 auto& pendant_edge = pendant_node.link().edge();
244 auto& proximal_edge = pendant_edge.primary_link().next().edge();
245 auto& distal_edge = pendant_edge.primary_link().next().next().edge();
248 assert(
degree( pendant_edge.primary_node() ) == 3 );
249 assert(
is_leaf( pendant_edge.secondary_node() ));
258 assert( pend_data.branch_length == 0.0 );
262 assert( placement_pair.placement->proximal_length >= used_length );
265 pend_data.branch_length = placement_pair.placement->pendant_length;
266 prox_data.
branch_length = placement_pair.placement->proximal_length - used_length;
270 node_data.
name = name_prefix + pquery_name.name;
273 used_length = placement_pair.placement->proximal_length;
274 insertion_edge = &distal_edge;
288 assert( end_edge_data.branch_length == 0.0 );
289 if( used_length < orig_length ) {
290 end_edge_data.branch_length = orig_length - used_length;
300 auto process_edge_multifurcating = [ &name_prefix ] (
303 std::vector<PlacementPair>
const& placement_pairs
309 auto& base_edge = new_node.link().edge();
310 auto& pri_edge = base_edge.primary_link().next().edge();
311 auto& sec_edge = base_edge.primary_link().next().next().edge();
314 assert(
degree( base_edge.primary_node() ) == 3 );
315 assert(
is_leaf( base_edge.secondary_node() ));
323 assert( base_data.branch_length == 0.0 );
327 auto const avg_prox_len = std::accumulate(
328 placement_pairs.begin(),
329 placement_pairs.end(),
331 [](
double in, PlacementPair
const& cur ){
332 return in + cur.placement->proximal_length;
334 ) /
static_cast<double>( placement_pairs.size() );
349 assert( placement_pairs.size() > 0 );
350 auto const min_pen_len = std::max( 0.0, std::min_element(
351 placement_pairs.begin(),
352 placement_pairs.end(),
353 []( PlacementPair
const& lhs, PlacementPair
const& rhs ){
354 return lhs.placement->pendant_length < rhs.placement->pendant_length;
356 )->placement->pendant_length );
360 assert( min_pen_len == placement_pairs[0].placement->pendant_length );
363 base_data.branch_length = min_pen_len;
366 for(
auto const& placement_pair : placement_pairs ) {
368 for(
auto const& pquery_name : placement_pair.pquery->names() ) {
372 auto& p_edge = p_node.link().edge();
377 assert( placement_pair.placement->pendant_length >= min_pen_len );
378 p_data.branch_length = placement_pair.placement->pendant_length - min_pen_len;
382 p_node_data.
name = name_prefix + pquery_name.name;
392 for(
size_t edge_i = 0; edge_i < edge_list.size(); ++edge_i ) {
394 assert( edge_list[ edge_i ] );
397 if( place_map[ edge_i ].size() == 0 ) {
403 if( place_map[ edge_i ].size() == 1 && place_map[ edge_i ][0].pquery->name_size() <= 1 ) {
404 assert( place_map[ edge_i ][0].pquery && place_map[ edge_i ][0].placement );
405 add_lonely_placement( result, *edge_list[ edge_i ], place_map[ edge_i ][0] );
410 if( fully_resolve ) {
411 process_edge_fully_resolved(
412 result, *edge_list[ edge_i ], place_map[ edge_i ]
415 process_edge_multifurcating(
416 result, *edge_list[ edge_i ], place_map[ edge_i ]