1 #include "../lib/bitvector.hpp"
2 #include "../lib/io_utils.hpp"
3 #include "../lib/multitree_iterator.hpp"
4 #include "../lib/subtree_extraction_impl.hpp"
5 #include "../lib/supertree_enumerator.hpp"
6 #include "../lib/supertree_variants.hpp"
7 #include "../lib/supertree_variants_multitree.hpp"
8 #include "../lib/validation.hpp"
9 #include <iostream>
10 #include <terraces/constraints.hpp>
11 #include <terraces/parser.hpp>
12 #include <terraces/rooting.hpp>
13
14 enum class mode { count, enumerate };
15
16 using namespace terraces;
17
is_comprehensive(const tree & t,const bitmatrix & matrix)18 bool is_comprehensive(const tree& t, const bitmatrix& matrix) {
19 const auto node_occ = compute_node_occ(t, matrix).first;
20 for (terraces::index i = 0; i < matrix.cols(); ++i) {
21 if (induced_lca(t, node_occ, i) != 0) {
22 return false;
23 }
24 }
25 return true;
26 }
27
count_and_check_trees(const multitree_node * multitree,const bitmatrix & matrix,const std::vector<tree> & ref_trees,const name_map & names)28 void count_and_check_trees(const multitree_node* multitree, const bitmatrix& matrix,
29 const std::vector<tree>& ref_trees, const name_map& names) {
30 std::cerr << "Enumerating all trees" << std::endl;
31
32 // store all supertrees
33 std::vector<tree> supertrees;
34 multitree_iterator mit(multitree);
35 do {
36 supertrees.push_back(mit.tree());
37 } while (mit.next());
38 if (big_integer{supertrees.size()} != multitree->num_trees) {
39 std::cerr << "PANIC! Inconsistent tree count" << std::endl;
40 }
41
42 // check for supertree property
43 for (const auto& supertree : supertrees) {
44 auto trees = subtrees(supertree, matrix);
45 for (terraces::index i = 0; i < trees.size(); ++i) {
46 if (!is_isomorphic_rooted(ref_trees[i], trees[i])) {
47 std::cerr << "PANIC! Invalid supertree:\nPartition " << i
48 << "\nSupertree: " << as_newick(supertree, names)
49 << "\nExpected: " << as_newick(ref_trees[i], names)
50 << "\nReceived: " << as_newick(trees[i], names)
51 << std::endl;
52 }
53 }
54 }
55
56 // check for isomorphic trees
57 std::vector<std::pair<std::vector<simple_bitvector>, terraces::index>> tree_bips;
58 for (terraces::index i = 0; i < supertrees.size(); ++i) {
59 const auto& supertree = supertrees[i];
60 tree_bips.emplace_back(tree_bipartitions(supertree), i);
61 }
62 std::sort(tree_bips.begin(), tree_bips.end());
63 int dupes = 0;
64 for (terraces::index i = 1; i < supertrees.size(); ++i) {
65 if (tree_bips[i - 1].first == tree_bips[i].first) {
66 ++dupes;
67 std::cerr << "Duplicate tree found:\n"
68 << as_newick(supertrees[tree_bips[i].second], names) << '\n'
69 << as_newick(supertrees[tree_bips[i - 1].second], names)
70 << std::endl;
71 }
72 }
73 std::cout << (supertrees.size() - dupes) << "(" << supertrees.size() << ")";
74 }
75
main(int argc,char ** argv)76 int main(int argc, char** argv) {
77 std::string usage = argc < 3 ? "" : argv[1];
78 mode m;
79 if (usage == "count") {
80 m = mode::count;
81 } else if (usage == "enumerate") {
82 m = mode::enumerate;
83 } else {
84 std::cerr << "Usage: " << argv[0] << " {count|enumerate} [treefile] [datafile]\n";
85 return -1;
86 }
87
88 try {
89 auto data_stream = utils::open_ifstream(argv[3]);
90 auto data = parse_bitmatrix(data_stream);
91 auto tree = parse_nwk(utils::read_file_full(argv[2]), data.indices);
92 auto num_leaves = data.matrix.rows();
93
94 for (terraces::index i = 1; i < tree.size(); ++i) {
95 std::cerr << "Rooting at node " << i << ":\n";
96 auto rerooted_tree = reroot_at_node(tree, i);
97 std::cerr << as_newick(rerooted_tree, data.names) << std::endl;
98 auto subtrees = terraces::subtrees(rerooted_tree, data.matrix);
99 auto constraints = compute_constraints(subtrees);
100 std::cerr << "Extracted subtrees and constraints\n";
101 std::cerr << "Removed " << deduplicate_constraints(constraints)
102 << " duplicate constraints" << std::endl;
103 auto root_split = terraces::root_split(rerooted_tree, num_leaves);
104 if (m == mode::count) {
105 tree_enumerator<variants::count_callback<big_integer>> e{{}};
106 auto c1 = e.run(num_leaves, constraints, root_split);
107 auto c2 = e.run(num_leaves, constraints);
108 std::cerr << "Counted supertrees" << std::endl;
109 std::cout << c1 << "\t" << c2 << std::endl;
110 } else {
111 tree_enumerator<variants::multitree_callback> e{{}};
112 auto t1 = e.run(num_leaves, constraints, root_split);
113 auto t2 = e.run(num_leaves, constraints);
114 std::cerr << "Checking supertrees" << std::endl;
115 std::cout << is_comprehensive(rerooted_tree, data.matrix) << "\t";
116 count_and_check_trees(t1, data.matrix, subtrees, data.names);
117 std::cout << "\t";
118 count_and_check_trees(t2, data.matrix, subtrees, data.names);
119 std::cout << std::endl;
120 }
121 }
122 } catch (std::exception& e) {
123 std::cerr << "Error: " << e.what() << "\n";
124 }
125 }