placement/simulator/functions.cpp

Go to the documentation of this file.

/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2019 Lucas Czech and HITS gGmbH
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
    Contact:
    Lucas Czech <lucas.czech@h-its.org>
    Exelixis Lab, Heidelberg Institute for Theoretical Studies
    Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
*/
 
#include "genesis/placement/simulator/functions.hpp"
 
#include "genesis/placement/function/functions.hpp"
#include "genesis/placement/function/helper.hpp"
#include "genesis/placement/function/masses.hpp"
 
#include "genesis/placement/simulator/distributions.hpp"
 
#include "genesis/tree/function/distances.hpp"
#include "genesis/tree/function/functions.hpp"
 
#include "genesis/utils/containers/matrix.hpp"
#include "genesis/utils/core/std.hpp"
#include "genesis/utils/text/string.hpp"
 
#include <cassert>
#include <cmath>
#include <ostream>
#include <stdexcept>
#include <vector>
 
namespace genesis {
namespace placement {
 
// =================================================================================================
//     Stream Output Operators
// =================================================================================================
 
std::ostream& operator << ( std::ostream& out, SimulatorEdgeDistribution const& distrib )
{
    out << "Weight of each edge: ";
    out << utils::join( distrib.edge_weights, ", ") << "\n";
    return out;
}
 
std::ostream& operator << ( std::ostream& out, SimulatorExtraPlacementDistribution const& distrib )
{
    out << "Extra placement weights:\n";
    for( size_t i = 0; i < distrib.placement_number_weights.size(); ++i ) {
        out << i << ": " << distrib.placement_number_weights[i] << "\n";
    }
    out << "Path length weights:\n";
    for( size_t i = 0; i < distrib.placement_path_length_weights.size(); ++i ) {
        out << i << ": " << distrib.placement_path_length_weights[i] << "\n";
    }
    return out;
}
 
std::ostream& operator << ( std::ostream& out, SimulatorLikeWeightRatioDistribution const& distrib )
{
    if( distrib.intervals.size() != distrib.weights.size() ) {
        throw std::logic_error( "Invalid SimulatorLikeWeightRatioDistribution." );
    }
 
    out << "Like weight ratio intervals and weights:\n";
    for( size_t i = 0; i < distrib.intervals.size(); ++i ) {
        out << distrib.intervals[i] << ": " << distrib.weights[i] << "\n";
    }
    return out;
}
 
// =================================================================================================
//     Edge Distribution : Set Weights
// =================================================================================================
 
// -----------------------------------------------------
//     set_uniform_weights
// -----------------------------------------------------
 
void set_uniform_weights( Sample const& sample, SimulatorEdgeDistribution& edge_distrib )
{
    set_uniform_weights( sample.tree().edge_count(), edge_distrib );
}
 
void set_uniform_weights( size_t edge_count, SimulatorEdgeDistribution& edge_distrib )
{
    edge_distrib.edge_weights = std::vector<double>( edge_count, 1.0 );
}
 
// -----------------------------------------------------
//     set_random_weights
// -----------------------------------------------------
 
void set_random_weights( Sample const& sample, SimulatorEdgeDistribution& edge_distrib )
{
    set_random_weights( sample.tree().edge_count(), edge_distrib );
}
 
void set_random_weights( size_t edge_count, SimulatorEdgeDistribution& edge_distrib )
{
    edge_distrib.edge_weights = std::vector<double>( edge_count, 0.0 );
 
    std::uniform_real_distribution<double> distrib( 0.0, 1.0 );
    for( size_t i = 0; i < edge_count; ++i ) {
        edge_distrib.edge_weights[i] = distrib( utils::Options::get().random_engine() );
    }
}
 
// -----------------------------------------------------
//     set_random_edges
// -----------------------------------------------------
 
void set_random_edges( Sample const& sample, SimulatorEdgeDistribution& edge_distrib )
{
    set_random_edges( sample.tree().edge_count(), edge_distrib );
}
 
void set_random_edges( size_t edge_count, SimulatorEdgeDistribution& edge_distrib )
{
    edge_distrib.edge_weights = std::vector<double>( edge_count, 0.0 );
 
    std::bernoulli_distribution distrib;
    for (size_t i = 0; i < edge_count; ++i) {
        if( distrib( utils::Options::get().random_engine() )) {
            edge_distrib.edge_weights[i] = 1.0;
        }
    }
}
 
// -----------------------------------------------------
//     set_depths_distributed_weights
// -----------------------------------------------------
 
void set_depths_distributed_weights( Sample const& sample, SimulatorEdgeDistribution& edge_distrib )
{
    auto depth_weights = closest_leaf_weight_distribution( sample );
    set_depths_distributed_weights( sample, depth_weights, edge_distrib );
}
 
void set_depths_distributed_weights(
    Sample const&              sample,
    std::vector<double> const& depth_weights,
    SimulatorEdgeDistribution& edge_distrib
) {
    // TODO Some of the code is copied from Sample::closest_leaf_weight_distribution(). Maybe
    //      it is worth putting this into a method which returns the closest leaf for edges instead
    //      of nodes.
 
    // Prepare weights vector.
    size_t num_edges = sample.tree().edge_count();
    edge_distrib.edge_weights = std::vector<double>(num_edges, 0.0);
 
    // Get a vector telling us the depth from each node to its closest leaf node.
    auto depths = closest_leaf_depth_vector(sample.tree());
 
    // Set the weight of each edge according to its depth in the tree.
    for( auto const& edge : sample.tree().edges() ) {
        // Try both nodes at the end of the edge and see which one is closer to a leaf.
        auto dp = depths[ edge.primary_node().index() ].second;
        auto ds = depths[ edge.secondary_node().index() ].second;
        auto ld = std::min(dp, ds);
 
        // Some safty. This holds as long as the indices are correct.
        assert( edge.index() < num_edges );
 
        // If the depth of the current edge is in the depth vector, use it.
        // Otherwise, the tree is deeper than the given depth vector, so use zero instead,
        // which will result in no placements being generated on this edge.
        if (ld < depth_weights.size()) {
            edge_distrib.edge_weights[ edge.index() ] = depth_weights[ld];
        } else {
            edge_distrib.edge_weights[ edge.index() ] = 0.0;
        }
    }
}
 
// -----------------------------------------------------
//     set_random_subtree_weights
// -----------------------------------------------------
 
size_t set_random_subtree_weights( Sample const& sample, SimulatorEdgeDistribution& edge_distrib )
{
    if( sample.tree().empty() ) {
        throw std::runtime_error("Cannot set weights using an empty Tree.");
    }
 
    // Reset all edge weights.
    size_t edge_count = sample.tree().edge_count();
    edge_distrib.edge_weights = std::vector<double>( edge_count, 0.0 );
 
    assert( edge_count > 0 );
    std::uniform_int_distribution<size_t> edge_selection( 0, edge_count - 1 );
    size_t edge_idx = edge_selection( utils::Options::get().random_engine() );
 
    // Randomly select a subtree, either starting at the end of the edge in direction of the root,
    // or away from it.
    PlacementTreeLink const* start_link;
    std::bernoulli_distribution dir_distrib;
    if( dir_distrib( utils::Options::get().random_engine() )) {
        // Primary direction
        start_link = &sample.tree().edge_at(edge_idx).primary_link();
    } else {
        // Secondary direction
        start_link = &sample.tree().edge_at(edge_idx).secondary_link();
    }
 
    for (
        auto cur_link = &start_link->next();
        cur_link != start_link;
        cur_link = &cur_link->outer().next()
    ) {
        edge_distrib.edge_weights[ cur_link->edge().index() ] = 1.0;
    }
 
    return edge_idx;
}
 
// -----------------------------------------------------
//     set_subtree_weights
// -----------------------------------------------------
 
void set_subtree_weights(
    Sample const& sample,
    size_t link_index,
    SimulatorEdgeDistribution& edge_distrib
) {
    // Validity checks.
    if( link_index >= sample.tree().link_count() ) {
        throw std::runtime_error( "Invalid link index for subtree." );
    }
    if( is_leaf( sample.tree().link_at( link_index ) )) {
        throw std::runtime_error( "Cannot use a leaf node as subtree." );
    }
 
    // Reset all edge weights.
    auto edge_count = sample.tree().edge_count();
    edge_distrib.edge_weights = std::vector<double>( edge_count, 0.0 );
 
    // Iterate the subtree and set edge weights.
    PlacementTreeLink const* start_link = &sample.tree().link_at( link_index );
    for (
        auto cur_link = &start_link->next();
        cur_link != start_link;
        cur_link = &cur_link->outer().next()
    ) {
        edge_distrib.edge_weights[ cur_link->edge().index() ] = 1.0;
    }
}
 
// -----------------------------------------------------
//     learn_per_edge_weights
// -----------------------------------------------------
 
void learn_per_edge_weights( Sample const& sample, SimulatorEdgeDistribution& edge_distrib )
{
    edge_distrib.edge_weights = placement_mass_per_edges_with_multiplicities( sample );
}
 
// =================================================================================================
//     Placement Number Distribution
// =================================================================================================
 
void learn_placement_number_weights(
    Sample const& sample,
    SimulatorExtraPlacementDistribution& p_distib
) {
    auto weights = std::vector<double>();
    for( auto const& pquery : sample.pqueries() ) {
        size_t extra_placements = pquery.placement_size() - 1;
        if( weights.size() < extra_placements + 1 ) {
            weights.resize( extra_placements + 1 );
        }
        ++weights[ extra_placements ];
    }
    p_distib.placement_number_weights = weights;
}
 
void learn_placement_path_length_weights(
    Sample const& sample,
    SimulatorExtraPlacementDistribution& p_distib
) {
    // The distance (path length) between two placements is the number of nodes between them.
    // Get a matrix that gives us this number for each pair of edges of the tree.
    auto edge_dist_matrix = tree::edge_path_length_matrix( sample.tree() );
 
    // Iterate all pqueries and collect the distances between all of their placements, where
    // the distance is the number of nodes between them.
    auto weights = std::vector<double>();
    for( auto const& pquery : sample.pqueries() ) {
        for( auto const& place_from : pquery.placements() ) {
            for( auto const& place_to : pquery.placements() ) {
                size_t dist = edge_dist_matrix(
                    place_from.edge().index(), place_to.edge().index()
                );
 
                // The following assertion holds as long as our dist matrix is correct.
                // If assertions are disabled (release), the whole condition should be eliminated
                // from the binary, so there is no speed penalty.
                if( place_from.edge().index() == place_to.edge().index() ) {
                    assert( dist == 0 );
                }
 
                // We don't need paths of length 0 currently. They will be eliminated in the
                // distrubution class anyway, when calling prepare().
                if( dist == 0 ) {
                    continue;
                }
                if( weights.size() < dist + 1 ) {
                    weights.resize( dist + 1 );
                }
                weights[ dist ] += 1.0;
            }
        }
    }
    p_distib.placement_path_length_weights = weights;
}
 
// =================================================================================================
//     Like Weight Ratio Distribution
// =================================================================================================
 
void learn_like_weight_ratio_distribution(
    Sample const& sample,
    SimulatorLikeWeightRatioDistribution& lwr_distib,
    size_t number_of_intervals
) {
    // Init the result vectors (we need one more element than number of intervals).
    auto intervals = std::vector<double>( number_of_intervals + 1, 0.0 );
    auto weights   = std::vector<double>( number_of_intervals + 1, 0.0 );
    for( size_t i = 0; i < number_of_intervals + 1; ++i ) {
        intervals[i] = static_cast<double>( i ) / static_cast<double>( number_of_intervals );
    }
 
    // Iterate all placements and use their rounded lwr value to increase the weights of the
    // distribution.
    for( auto const& pquery : sample.pqueries() ) {
        for( auto const& placement : pquery.placements() ) {
 
            double rounded_lwr = std::round( placement.like_weight_ratio * number_of_intervals )
                               / number_of_intervals;
 
            if( 0.0 > rounded_lwr || rounded_lwr > 1.0 ) {
                throw std::runtime_error( "Invalid like_weight_ratio in Sample." );
            }
            size_t pos = static_cast<size_t>( rounded_lwr * number_of_intervals );
            assert( pos < weights.size() );
            weights[ pos ] += 1.0;
        }
    }
 
    // Set the result.
    lwr_distib.intervals = intervals;
    lwr_distib.weights   = weights;
}
 
} // namespace placement
} // namespace genesis