63 double const scaler_lhs,
64 double const scaler_rhs
74 double const scaler_lhs,
75 double const scaler_rhs
79 throw std::runtime_error(
"Incompatible MassTrees for merging." );
83 if( scaler_lhs != 1.0 ) {
87 #pragma omp parallel for 88 for(
size_t i = 0; i < lhs.
edge_count(); ++i ) {
91 lhs_masses[ rhs_mass.first ] += scaler_rhs * rhs_mass.second;
98 #pragma omp parallel for 99 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
106 #pragma omp parallel for 107 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
109 for(
auto& mass : masses ) {
110 mass.second = -mass.second;
117 #pragma omp parallel for 118 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
120 for(
auto& mass : masses ) {
121 mass.second *= factor;
130 if( total_mass == 0.0 ) {
134 #pragma omp parallel for 135 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
137 mass.second /= total_mass;
144 #pragma omp parallel for 145 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
147 std::map<double, double> relative;
149 for(
auto& mass : edge_data.masses ) {
150 relative[ mass.first / edge_data.branch_length ] += mass.second;
153 edge_data.
masses = relative;
154 edge_data.branch_length = 1.0;
162 #pragma omp parallel for 163 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
167 double central_mass = 0.0;
169 for(
auto const& mass : edge_data.masses ) {
171 work += mass.second * std::abs( branch_center - mass.first );
172 central_mass += mass.second;
175 edge_data.masses.clear();
176 edge_data.masses[ branch_center ] = central_mass;
185 #pragma omp parallel for 186 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
192 if( edge_data.masses.empty() ) {
196 double mass_center = 0.0;
197 double mass_sum = 0.0;
201 for(
auto const& mass : edge_data.masses ) {
202 mass_center += mass.first * mass.second;
203 mass_sum += mass.second;
207 mass_center /= mass_sum;
210 for(
auto const& mass : edge_data.masses ) {
212 work += mass.second * std::abs( mass_center - mass.first );
216 edge_data.masses.clear();
217 edge_data.masses[ mass_center ] = mass_sum;
224 if( number_of_bins == 0 ) {
225 throw std::invalid_argument(
"Cannot use number_of_bins == 0." );
231 auto get_bin_pos = [ number_of_bins ](
double pos,
double bl )
233 auto const nb =
static_cast<double>( number_of_bins );
236 auto const pn = std::min( std::max( 0.0, pos / bl * nb ), std::nextafter( nb, 0.0 ));
240 return ( std::floor( pn ) * bl / nb ) + ( bl / nb / 2.0 );
245 #pragma omp parallel for 246 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
250 auto new_masses = std::map<double, double>();
253 for(
auto const& mass : edge_data.masses ) {
254 auto const bin = get_bin_pos( mass.first, edge_data.branch_length );
256 work += mass.second * std::abs( bin - mass.first );
257 new_masses[ bin ] += mass.second;
261 edge_data.
masses = new_masses;
274 if( mass_trees.size() < 2 ) {
279 size_t const num_edges = mass_trees[0].edge_count();
280 auto avg_br_lens = std::vector<double>( num_edges, 0.0 );
283 for(
auto const& tree : mass_trees ) {
286 if( tree.edge_count() != num_edges ) {
287 throw std::runtime_error(
288 "Cannot make average branch lengths, because trees have different size." 293 assert( avg_br_lens.size() == tree.edge_count() );
294 for(
size_t edge_idx = 0; edge_idx < num_edges; ++edge_idx ) {
295 avg_br_lens[edge_idx] += tree.edge_at( edge_idx ).data<
MassTreeEdgeData>().branch_length;
300 for(
auto& ae : avg_br_lens ) {
301 ae /=
static_cast<double>( mass_trees.size() );
305 for(
auto& tree : mass_trees ) {
307 #pragma omp parallel for 308 for(
size_t edge_idx = 0; edge_idx < tree.edge_count(); ++edge_idx ) {
312 auto new_masses = std::map<double, double>();
315 auto const scaler = avg_br_lens[ edge_idx ] / edge_data.
branch_length;
318 for(
auto const& mass : edge_data.masses ) {
320 auto const new_pos = mass.first * scaler;
321 new_masses[ new_pos ] += mass.second;
325 edge_data.
masses = new_masses;
333 auto result = std::vector<double>( tree.
edge_count(), 0.0 );
335 #pragma omp parallel for 336 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
337 auto const& edge = tree.
edge_at(i);
338 assert( i == edge.index() );
352 if( mass_trees.empty() || mass_trees[0].empty() ) {
358 #pragma omp parallel for 359 for(
size_t i = 0; i < mass_trees.size(); ++i ) {
360 if( mass_trees[i].
edge_count() != result.cols() ) {
361 throw std::runtime_error(
362 "Cannot calculate masses per edge for a Tree set with Trees " 363 "with unequal edge count." 376 auto result = std::vector<std::pair<double, double>>( tree.
edge_count(), { 0.0, 0.0 });
378 #pragma omp parallel for 379 for(
size_t i = 0; i < tree.
edge_count(); ++i ) {
380 auto const& edge = tree.
edge_at(i);
384 if( edge_data.masses.empty() ) {
389 double mass_pos = 0.0;
390 double mass_sum = 0.0;
391 for(
auto const& mass : edge_data.masses ) {
392 mass_pos += mass.first * mass.second;
393 mass_sum += mass.second;
397 mass_pos /= mass_sum;
399 result[ edge.index() ].first = mass_pos;
400 result[ edge.index() ].second = mass_sum;
408 double total_mass = 0.0;
409 for(
auto const& edge : tree.
edges() ) {
411 total_mass += mass.second;
421 LOG_INFO <<
"Invalid Tree topology.";
424 if( ! tree_data_is< MassTreeNodeData, MassTreeEdgeData >( tree )) {
425 LOG_INFO <<
"Tree does not only contain Mass Node and Edge data types.";
430 double mass_sum = 0.0;
431 for(
auto const& edge : tree.
edges() ) {
432 auto const edge_data =
dynamic_cast< MassTreeEdgeData const*
>( edge.data_ptr() );
433 if( edge_data ==
nullptr ) {
434 LOG_INFO <<
"Edge data type is not 'MassTreeEdgeData'.";
438 for(
auto const& mass : edge_data->masses ) {
439 if( mass.first < 0.0 ) {
440 LOG_INFO <<
"Mass with branch position < 0.0";
443 if( mass.first > edge_data->branch_length ) {
444 LOG_INFO <<
"Mass with branch position > branch_length";
448 mass_sum += mass.second;
452 if( valid_total_mass_difference >= 0.0 && std::abs(mass_sum) > valid_total_mass_difference ) {
453 LOG_INFO <<
"Total mass difference " << mass_sum
454 <<
" is higher than " << valid_total_mass_difference;
void mass_tree_reverse_signs(MassTree &tree)
Reverse the sign of each mass point on a MassTree.
Data class for MassTreeEdges. Stores the branch length and a list of masses with their positions alon...
double mass_tree_center_masses_on_branches_averaged(MassTree &tree)
Accumulate all masses of a MassTree at the average mass position per edge.
size_t edge_count() const
Return the number of TreeEdges of the Tree.
double sum(const Histogram &h)
double mass_tree_sum_of_masses(MassTree const &tree)
Return the total sum of all masses on the MassTree.
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
bool mass_tree_validate(MassTree const &tree, double valid_total_mass_difference)
Validate the data on a MassTree.
std::vector< double > mass_tree_mass_per_edge(MassTree const &tree)
Return a std::vector that contains the total mass for each edge of the given MassTree.
std::map< double, double > masses
List of masses stored on this branch, sorted by their position on the branch.
size_t edge_count(Tree const &tree)
Return the number of Edges of a Tree. Same as Tree::edge_count().
void mass_tree_transform_to_unit_branch_lengths(MassTree &tree)
Set all branch lengths of a MassTree to 1.0, while keeping the relative position of all masses on the...
std::vector< std::pair< double, double > > mass_tree_mass_per_edge_averaged(MassTree const &tree)
Return a std::vector that contains the total Mass for each edge of the given MassTree (second value)...
Class for representing phylogenetic trees.
void mass_tree_merge_trees_inplace(MassTree &lhs, MassTree const &rhs, double const scaler_lhs, double const scaler_rhs)
Merge all masses of two MassTrees by adding them to the first MassTree.
utils::Range< IteratorEdges > edges()
void mass_trees_make_average_branch_lengths(std::vector< MassTree > &mass_trees)
Change the branch lengths of all trees to their average, and move the masses accordingly in a proport...
TreeEdge & edge_at(size_t index)
Return the TreeEdge at a certain index.
Provides easy and fast logging functionality.
void mass_tree_clear_masses(MassTree &tree)
Clear all masses of a MassTree, while keeping its topology.
double branch_length
Branch length of the edge.
bool validate_topology(Tree const &tree)
Validate that all internal pointers of the Tree elements (TreeLinks, TreeNodes, TreeEdges) to each ot...
void mass_tree_scale_masses(MassTree &tree, double factor)
Scale all masses of a MassTree with the multiplicative factor factor.
void mass_tree_normalize_masses(MassTree &tree)
Scale all masses of a MassTree so that they sum up to 1.0.
MassTree mass_tree_merge_trees(MassTree const &lhs, MassTree const &rhs, double const scaler_lhs, double const scaler_rhs)
Merge all masses of two MassTrees into one and return it.
double mass_tree_center_masses_on_branches(MassTree &tree)
Accumulate all masses of a MassTree on the centers of their edges.
double mass_tree_binify_masses(MassTree &tree, size_t number_of_bins)
Accumulate all masses of a MassTree into bins on each branch.
#define LOG_INFO
Log an info message. See genesis::utils::LoggingLevel.