common.hpp

Go to the documentation of this file.

#ifndef GENESIS_UTILS_MATH_COMMON_H_
#define GENESIS_UTILS_MATH_COMMON_H_
 
/*
    Genesis - A toolkit for working with phylogenetic data.
    Copyright (C) 2014-2021 Lucas Czech
 
    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 <lczech@carnegiescience.edu>
    Department of Plant Biology, Carnegie Institution For Science
    260 Panama Street, Stanford, CA 94305, USA
*/
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <functional>
#include <limits>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace genesis {
namespace utils {
 
// =================================================================================================
//     Constants and General Functions
// =================================================================================================
 
constexpr double PI = 3.141592653589793238463;
 
inline double circumference( double radius )
{
    return 2.0 * PI * radius;
}
 
double log_factorial( size_t n );
 
size_t binomial_coefficient( size_t n, size_t k );
 
double binomial_coefficient_approx( size_t n, size_t k, bool lenient = false );
 
inline double binomial_distribution( size_t k, size_t n, double p, bool lenient = false )
{
    if( ! std::isfinite(p) || p < 0.0 || p > 1.0 ) {
        throw std::invalid_argument(
            "Cannot compute binomial distribution with p outside of [ 0, 1 ]"
        );
    }
    double const coeff = binomial_coefficient_approx( n, k, lenient );
    return coeff * std::pow( p, k ) * std::pow( 1.0 - p, n - k );
}
 
// =================================================================================================
//     Number Handling
// =================================================================================================
 
template< typename T > inline constexpr
T abs_diff( T const& lhs, T const& rhs )
{
    return ((lhs > rhs) ? (lhs - rhs) : (rhs - lhs));
}
 
template <typename T> inline constexpr
int signum(T x, std::false_type )
{
    // The tag type `std::false_type` is unnamed in order to avoid compiler warnings about an
    // unused parameter. As this function is `constexpr`, we cannot shut this warning down by
    // using `(void) param`, so making it unnamed is a reasonable solution in this case.
    return T(0) < x;
}
 
template <typename T> inline constexpr
int signum(T x, std::true_type )
{
    // The tag type `std::false_type` is unnamed in order to avoid compiler warnings about an
    // unused parameter. As this function is `constexpr`, we cannot shut this warning down by
    // using `(void) param`, so making it unnamed is a reasonable solution in this case.
    return (T(0) < x) - (x < T(0));
}
 
template <typename T> inline constexpr
int signum(T x)
{
    return signum( x, std::is_signed<T>() );
}
 
inline bool almost_equal_relative(
    double lhs,
    double rhs,
    double max_rel_diff = std::numeric_limits<double>::epsilon()
) {
    // Calculate the difference.
    auto diff = std::abs( lhs - rhs );
 
    // Find the larger number.
    auto largest = std::max( std::abs( lhs ), std::abs( rhs ));
 
    // Do the comparison.
    return ( diff <= largest * max_rel_diff );
}
 
inline double round_to( double x, size_t accuracy_order )
{
    double factor = std::pow( 10, accuracy_order );
    return std::round( x * factor ) / factor;
}
 
inline size_t int_pow( size_t base, size_t exp )
{
    // Using Exponentiation by squaring, see
    // http://stackoverflow.com/a/101613/4184258
    size_t result = 1;
    while( exp )  {
        if( exp & 1 ) {
            result *= base;
        }
        exp >>= 1;
        base *= base;
    }
    return result;
}
 
inline bool is_valid_int_pow( size_t base, size_t exp )
{
    return std::pow( base, exp ) < static_cast<double>( std::numeric_limits<size_t>::max() );
}
 
inline constexpr double squared( double x )
{
    return x * x;
}
 
inline constexpr double cubed( double x )
{
    return x * x * x;
}
 
// =================================================================================================
//     Helper Functions
// =================================================================================================
 
template <class ForwardIteratorA, class ForwardIteratorB>
std::pair<std::vector<double>, std::vector<double>> finite_pairs(
    ForwardIteratorA first_a, ForwardIteratorA last_a,
    ForwardIteratorB first_b, ForwardIteratorB last_b
) {
    // Prepare result.
    std::vector<double> vec_a;
    std::vector<double> vec_b;
 
    // Iterate in parallel.
    while( first_a != last_a && first_b != last_b ) {
        if( std::isfinite( *first_a ) && std::isfinite( *first_b ) ) {
            vec_a.push_back( *first_a );
            vec_b.push_back( *first_b );
        }
        ++first_a;
        ++first_b;
    }
    if( first_a != last_a || first_b != last_b ) {
        throw std::runtime_error(
            "Ranges need to have same length."
        );
    }
 
    assert( vec_a.size() == vec_b.size() );
    return { vec_a, vec_b };
}
 
template <class ForwardIteratorA, class ForwardIteratorB>
void for_each_finite_pair(
    ForwardIteratorA first_a, ForwardIteratorA last_a,
    ForwardIteratorB first_b, ForwardIteratorB last_b,
    std::function<void( double value_a, double value_b )> execute
) {
    // Iterate in parallel.
    while( first_a != last_a && first_b != last_b ) {
        if( std::isfinite( *first_a ) && std::isfinite( *first_b ) ) {
            execute( *first_a, *first_b );
        }
        ++first_a;
        ++first_b;
    }
 
    // Check if the ranges have the same length.
    if( first_a != last_a || first_b != last_b ) {
        throw std::runtime_error( "Ranges need to have same length." );
    }
}
 
} // namespace utils
} // namespace genesis
 
#endif // include guard