1 #include <terraces/advanced.hpp>
2
3 #include <terraces/clamped_uint.hpp>
4 #include <terraces/errors.hpp>
5 #include <terraces/rooting.hpp>
6 #include <terraces/subtree_extraction.hpp>
7
8 #include "multitree_iterator.hpp"
9 #include "supertree_enumerator.hpp"
10 #include "supertree_variants.hpp"
11 #include "supertree_variants_multitree.hpp"
12
13 namespace terraces {
14
create_supertree_data(const tree & tree,const bitmatrix & data)15 supertree_data create_supertree_data(const tree& tree, const bitmatrix& data) {
16 auto root = find_comprehensive_taxon(data);
17 utils::ensure<bad_input_error>(data.rows() == num_leaves_from_nodes(tree.size()),
18 bad_input_error_type::tree_mismatching_size);
19 utils::ensure<no_usable_root_error>(root != none, "No comprehensive taxon found");
20 auto rerooted_tree = tree;
21 reroot_at_taxon_inplace(rerooted_tree, root);
22 auto trees = subtrees(rerooted_tree, data);
23 auto constraints = compute_constraints(trees);
24 deduplicate_constraints(constraints);
25
26 auto num_leaves = data.rows();
27 utils::ensure<bad_input_error>(num_leaves >= 4, bad_input_error_type::nwk_tree_trivial);
28 return {constraints, num_leaves, root};
29 }
30
find_comprehensive_taxon(const bitmatrix & data)31 index find_comprehensive_taxon(const bitmatrix& data) {
32 for (index i = 0; i < data.rows(); ++i) {
33 bool comp = true;
34 for (index j = 0; j < data.cols(); ++j) {
35 comp &= data.get(i, j);
36 }
37 if (comp) {
38 return i;
39 }
40 }
41 return none;
42 }
43
maximum_comprehensive_columnset(const bitmatrix & data)44 bitmatrix maximum_comprehensive_columnset(const bitmatrix& data) {
45 std::vector<index> row_counts(data.rows(), 0u);
46 for (index i = 0; i < data.rows(); ++i) {
47 for (index j = 0; j < data.cols(); ++j) {
48 row_counts[i] += data.get(i, j) ? 1u : 0u;
49 }
50 }
51 auto it = std::max_element(row_counts.begin(), row_counts.end());
52 index comp_row = static_cast<index>(std::distance(row_counts.begin(), it));
53 std::vector<index> columns;
54 for (index j = 0; j < data.cols(); ++j) {
55 if (data.get(comp_row, j)) {
56 columns.push_back(j);
57 }
58 }
59 return data.get_cols(columns);
60 }
61
fast_count_terrace(const supertree_data & data)62 index fast_count_terrace(const supertree_data& data) {
63 tree_enumerator<variants::check_callback> enumerator{{}};
64 try {
65 return enumerator.run(data.num_leaves, data.constraints, data.root);
66 } catch (terraces::tree_count_overflow_error&) {
67 return std::numeric_limits<index>::max();
68 }
69 }
70
check_terrace(const supertree_data & data)71 bool check_terrace(const supertree_data& data) { return fast_count_terrace(data) > 1; }
72
count_terrace(const supertree_data & data)73 index count_terrace(const supertree_data& data) {
74 tree_enumerator<variants::clamped_count_callback> enumerator{{}};
75 try {
76 return enumerator.run(data.num_leaves, data.constraints, data.root).value();
77 } catch (terraces::tree_count_overflow_error&) {
78 return std::numeric_limits<index>::max();
79 }
80 }
81
count_terrace_bigint(const supertree_data & data)82 big_integer count_terrace_bigint(const supertree_data& data) {
83 tree_enumerator<variants::count_callback<big_integer>> enumerator{{}};
84 return enumerator.run(data.num_leaves, data.constraints, data.root);
85 }
86
print_terrace_compressed(const supertree_data & data,const name_map & names,std::ostream & output)87 big_integer print_terrace_compressed(const supertree_data& data, const name_map& names,
88 std::ostream& output) {
89 tree_enumerator<variants::multitree_callback> enumerator{{}};
90 auto result = enumerator.run(data.num_leaves, data.constraints, data.root);
91 output << as_newick(result, names);
92
93 return result->num_trees;
94 }
95
print_terrace(const supertree_data & data,const name_map & names,std::ostream & output)96 big_integer print_terrace(const supertree_data& data, const name_map& names, std::ostream& output) {
97 tree_enumerator<variants::multitree_callback> enumerator{{}};
98 auto result = enumerator.run(data.num_leaves, data.constraints, data.root);
99 multitree_iterator mit{result};
100 do {
101 output << as_newick(mit.tree(), names) << '\n';
102 } while (mit.next());
103
104 return result->num_trees;
105 }
106
enumerate_terrace(const supertree_data & data,std::function<void (const tree &)> callback)107 void enumerate_terrace(const supertree_data& data, std::function<void(const tree&)> callback) {
108 tree_enumerator<variants::multitree_callback> enumerator{{}};
109 auto result = enumerator.run(data.num_leaves, data.constraints, data.root);
110 multitree_iterator mit{result};
111 do {
112 callback(mit.tree());
113 } while (mit.next());
114 }
115
116 } // namespace terraces
117