1 #ifndef SUPERTREE_VARIANTS_HPP
2 #define SUPERTREE_VARIANTS_HPP
3 
4 #include <terraces/constraints.hpp>
5 
6 #include <terraces/bigint.hpp>
7 #include <terraces/clamped_uint.hpp>
8 
9 #include "bipartitions.hpp"
10 #include "trees_impl.hpp"
11 
12 namespace terraces {
13 namespace variants {
14 
15 template <typename Result>
16 /**
17  * An abstract callback used to control the execution of \ref tree_enumerator.
18  * Every callback must declare a return type whose values are aggregated
19  * (sum over different bipartitions) and combined (results from left and right subtree).
20  * <p>Additionally, the enumeration can be terminated early using \ref fast_return
21  * or \ref continue_iteration, if one does not want to iterate over the whole recursion tree.</p>
22  */
23 class abstract_callback {
24 public:
25 	/** The result type returned by the @ref tree_enumerator. */
26 	using result_type = Result;
27 
28 	/**
29 	 * Called when a (sub)call begins.
30 	 * \param leaves The leaf set inspected at the current recursive call.
31 	 */
enter(const ranked_bitvector & leaves)32 	void enter(const ranked_bitvector& leaves) { (void)leaves; }
33 	/**
34 	 * Called when a (sub)call finishes.
35 	 * \param val The result that was either returned from \ref fast_return_value
36 	 *            or accumulated from the recursive subcalls.
37 	 * \returns The result that should be passed to the calling function.
38 	 */
exit(Result val)39 	Result exit(Result val) { return val; }
40 
41 	/** Returns the result for a single leaf. */
42 	Result base_one_leaf(index);
43 	/** Returns the result for two leaves. */
44 	Result base_two_leaves(index, index);
45 	/** Returns the result for multiple leaves without constraints. */
46 	Result base_unconstrained(const ranked_bitvector&);
47 	/** Returns an empty result. */
48 	Result null_result() const;
49 
50 	/**
51 	 * Called before iterating over the bipartitions to check if we can return early.
52 	 * This may be used to avoid searching subtrees that are 'not interesting'.
53 	 * \see fast_return_value
54 	 * \returns true if and only if we want to skip iterating over the current bipartitions.
55 	 *          (default: false)
56 	 */
fast_return(const bipartitions &)57 	bool fast_return(const bipartitions&) { return false; }
58 	/** Returns the result to be returned in case \ref fast_return is true. */
59 	Result fast_return_value(const bipartitions&);
60 
61 	/**
62 	 * Called when we begin iterating over the bipartitions.
63 	 * \param
64 	 * \returns The initial value for the result used for accumulation
65 	 *          - something equivalent to 0.
66 	 *          (default: value-initialization of @p Result)
67 	 */
begin_iteration(const bipartitions & bip_it,const bitvector & c_occ,const constraints & c)68 	Result begin_iteration(const bipartitions& bip_it, const bitvector& c_occ,
69 	                       const constraints& c) {
70 		(void)bip_it;
71 		(void)c_occ;
72 		(void)c;
73 		return Result{};
74 	}
75 	/**
76 	 * Returns true iff the iteration should continue.
77 	 * This may be used to stop iterating after a sufficient number of results have been
78 	 * accumulated.
79 	 * \returns true if and only if we want to continue iterating over
80 	 */
continue_iteration(Result)81 	bool continue_iteration(Result) { return true; }
82 	/** Called when an iteration step begins. */
step_iteration(const bipartitions &,index)83 	void step_iteration(const bipartitions&, index) {}
84 	/** Called when the last iteration step has finished. */
finish_iteration()85 	void finish_iteration() {}
86 	/** Called before descending into the left subset. */
left_subcall()87 	void left_subcall() {}
88 	/** Called before descending into the right subset. */
right_subcall()89 	void right_subcall() {}
90 	/**
91 	 * Accumulates the results from multiple bipartitions.
92 	 * \param accumulator The current value for the accumulator.
93 	 * \param value The value to be added to the accumulator.
94 	 * \returns The new accumulator value.
95 	 */
96 	Result accumulate(Result accumulator, Result value);
97 	/**
98 	 * Combines the results from two subcalls.
99 	 * \param left The result from the left subcall.
100 	 * \param right The result from the right subcall.
101 	 * \returns The combined result.
102 	 */
103 	Result combine(Result left, Result right);
104 };
105 
106 template <typename Number>
107 /**
108  * A callback implementation that counts all trees using the \p Number type.
109  * It gives correct results if the number type is integer and doesn't overflow.
110  */
111 class count_callback : public abstract_callback<Number> {
112 public:
113 	using return_type = typename abstract_callback<Number>::result_type;
114 	// only one choice for a single leaf
base_one_leaf(index)115 	return_type base_one_leaf(index) { return 1; }
116 	// only one choice for two leaves
base_two_leaves(index,index)117 	return_type base_two_leaves(index, index) { return 1; }
118 	// (#unrooted trees) choices for leaves with no constraints
base_unconstrained(const ranked_bitvector & leaves)119 	return_type base_unconstrained(const ranked_bitvector& leaves) {
120 		return count_unrooted_trees<return_type>(leaves.count());
121 	}
null_result() const122 	return_type null_result() const { return 0; }
123 
124 	// The number of bipartitions gives a lower bound on the number of trees.
fast_return_value(const bipartitions & bip_it)125 	index fast_return_value(const bipartitions& bip_it) { return bip_it.num_bip(); }
126 
127 	// Multiple choices are counted independently
accumulate(return_type acc,return_type val)128 	return_type accumulate(return_type acc, return_type val) { return acc + val; }
129 	// Choices from two the subtrees can be combined in any possible way
combine(return_type left,return_type right)130 	return_type combine(return_type left, return_type right) { return left * right; }
131 };
132 
133 /**
134  * A callback implementation that counts all trees using \ref clamped_uint.
135  * Thus, the result will be clamped at the maximal value for uint64_t,
136  * eliminating the need to fully walk through the recursion.
137  */
138 class clamped_count_callback : public count_callback<clamped_uint> {
139 public:
140 	using return_type = typename count_callback<clamped_uint>::result_type;
141 	// No need to keep counting if we already overflowed.
continue_iteration(return_type result)142 	bool continue_iteration(return_type result) { return !result.is_clamped(); }
143 };
144 
145 /**
146  * A callback implementation that returns a simple lower bound to the number
147  * of trees compatible with the given constraints.
148  * This allows to quickly check whether a tree lies on a phylogenetic terrace.
149  * It will always run faster than callbacks that enumerate all possible trees.
150  */
151 class check_callback : public abstract_callback<index> {
152 public:
153 	using return_type = index;
154 
155 	// No choices for a single leaf
base_one_leaf(index)156 	index base_one_leaf(index) { return 1; }
157 	// No choices for two leaves
base_two_leaves(index,index)158 	index base_two_leaves(index, index) { return 1; }
159 	// There are at least two possible trees if we have at least three unconstrained leaves
base_unconstrained(const ranked_bitvector &)160 	index base_unconstrained(const ranked_bitvector&) { return 2; }
null_result() const161 	return_type null_result() const { return 0; }
162 
163 	/* Since we assume there are no incompatible constraints,
164 	 * every subcall will return at least 1, so the number of bipartitions
165 	 * gives a simple lower bound for the number of trees. */
fast_return(const bipartitions & bip_it)166 	bool fast_return(const bipartitions& bip_it) { return bip_it.num_bip() > 1; }
fast_return_value(const bipartitions & bip_it)167 	index fast_return_value(const bipartitions& bip_it) { return bip_it.num_bip(); }
168 
169 	/* No need to keep on counting if we already know we are on a terrace. */
continue_iteration(index acc)170 	bool continue_iteration(index acc) { return acc < 2; }
171 
accumulate(index acc,index val)172 	index accumulate(index acc, index val) { return acc + val; }
combine(index left,index right)173 	index combine(index left, index right) { return left * right; }
174 };
175 
176 } // namespace variants
177 } // namespace terraces
178 
179 #endif // SUPERTREE_VARIANTS_HPP
180