distance.hpp

Go to the documentation of this file.

#ifndef GENESIS_UTILS_MATH_DISTANCE_H_
#define GENESIS_UTILS_MATH_DISTANCE_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 <lucas.czech@h-its.org>
    Exelixis Lab, Heidelberg Institute for Theoretical Studies
    Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
*/
 
#include "genesis/utils/core/algorithm.hpp"
#include "genesis/utils/containers/matrix.hpp"
#include "genesis/utils/math/common.hpp"
#include "genesis/utils/math/ranking.hpp"
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <stdexcept>
#include <utility>
#include <vector>
 
namespace genesis {
namespace utils {
 
// =================================================================================================
//     Norms
// =================================================================================================
 
template <class ForwardIterator>
double p_norm( ForwardIterator first, ForwardIterator last, double p = 2.0 )
{
    // Validity. We allow positive inifity.
    if( p < 1.0 || ( ! std::isfinite( p ) && ! std::isinf( p ))) {
        throw std::runtime_error( "Cannot calculate p-norm with p < 1.0" );
    }
    assert( p >= 1.0 );
    assert( std::isfinite( p ) || std::isinf( p ));
 
    double sum = 0.0;
    size_t cnt = 0;
 
    // Add vector elements.
    auto it = first;
    while( it != last ) {
        if( std::isfinite( *it ) ) {
            if( std::isfinite( p )) {
                sum += std::pow( std::abs( *it ), p );
            } else {
                sum = std::max( sum, std::abs( *it ));
            }
            ++cnt;
        }
        ++it;
    }
 
    // If there are no valid elements, return an all-zero result.
    if( cnt == 0 ) {
        return 0.0;
    }
 
    // Return the result.
    assert( cnt > 0 );
    if( std::isfinite( p )) {
        return std::pow( sum, 1.0 / p );
    } else {
        return sum;
    }
 
    // Make old compilers happy.
    return 0.0;
}
 
inline double p_norm( std::vector<double> const& vec, double p = 2.0 )
{
    return p_norm( vec.begin(), vec.end(), p );
}
 
template <class ForwardIterator>
double manhattan_norm( ForwardIterator first, ForwardIterator last )
{
    return p_norm( first, last, 1.0 );
}
 
inline double manhattan_norm( std::vector<double> const& vec )
{
    return p_norm( vec.begin(), vec.end(), 1.0 );
}
 
template <class ForwardIterator>
double euclidean_norm( ForwardIterator first, ForwardIterator last )
{
    return p_norm( first, last, 2.0 );
}
 
inline double euclidean_norm( std::vector<double> const& vec )
{
    return p_norm( vec.begin(), vec.end(), 2.0 );
}
 
template <class ForwardIterator>
double maximum_norm( ForwardIterator first, ForwardIterator last )
{
    return p_norm( first, last, std::numeric_limits<double>::infinity() );
}
 
inline double maximum_norm( std::vector<double> const& vec )
{
    return p_norm( vec.begin(), vec.end(), std::numeric_limits<double>::infinity() );
}
 
template <class ForwardIterator>
double aitchison_norm( ForwardIterator first, ForwardIterator last )
{
    double sum = 0.0;
    size_t cnt = 0;
 
    // Outer loop.
    auto it_out = first;
    while( it_out != last ) {
        if( std::isfinite( *it_out ) ) {
 
            if( *it_out <= 0.0 ) {
                throw std::invalid_argument(
                    "Cannot calculate Aitchison norm of non-positive values."
                );
            }
 
            // Inner loop.
            auto it_in = first;
            while( it_in != last ) {
                if( std::isfinite( *it_in ) ) {
                    auto const ln = std::log( *it_out / *it_in );
                    sum += ln * ln;
                }
                ++it_in;
            }
 
            ++cnt;
        }
        ++it_out;
    }
 
    // If there are no valid elements, return an all-zero result.
    if( cnt == 0 ) {
        return 0.0;
    }
 
    // Return the result.
    assert( cnt > 0 );
    return std::sqrt( sum / ( 2.0 * static_cast<double>( cnt )));
}
 
inline double aitchison_norm( std::vector<double> const& vec )
{
    return aitchison_norm( vec.begin(), vec.end() );
}
 
// =================================================================================================
//     Distances
// =================================================================================================
 
template <class ForwardIteratorA, class ForwardIteratorB>
double p_norm_distance(
    ForwardIteratorA first_a, ForwardIteratorA last_a,
    ForwardIteratorB first_b, ForwardIteratorB last_b,
    double p = 2.0
) {
    // Validity. We allow positive inifity.
    if( p < 1.0 || ( ! std::isfinite( p ) && ! std::isinf( p ))) {
        throw std::runtime_error( "Cannot calculate p-norm distance with p < 1.0" );
    }
    assert( p >= 1.0 );
    assert( std::isfinite( p ) || std::isinf( p ));
 
    // For "normal" p norms, just add up. For maximum norm (p=inf), we need a special case,
    // as double-precision arithmetics does not work the same as actual math :-)
    double sum = 0.0;
    size_t cnt = 0;
    if( std::isfinite( p )) {
        assert( p >= 1.0 );
        for_each_finite_pair( first_a, last_a, first_b, last_b, [&]( double val_a, double val_b ){
            sum += std::pow( std::abs( val_a - val_b ), p );
            ++cnt;
        });
    } else {
        assert( std::isinf( p ));
        for_each_finite_pair( first_a, last_a, first_b, last_b, [&]( double val_a, double val_b ){
            sum = std::max( sum, std::abs( val_a - val_b ) );
            ++cnt;
        });
    }
 
    // If there are no valid elements, return an all-zero result.
    if( cnt == 0 ) {
        return 0.0;
    }
 
    // Return the result.
    assert( cnt > 0 );
    if( std::isfinite( p )) {
        return std::pow( sum, 1.0 / p );
    } else {
        return sum;
    }
 
    // Make old compilers happy.
    return 0.0;
}
 
inline double p_norm_distance(
    std::vector<double> const& vec_a, std::vector<double> const& vec_b, double p = 2.0
) {
    return p_norm_distance( vec_a.begin(), vec_a.end(), vec_b.begin(), vec_b.end(), p );
}
 
template <class ForwardIteratorA, class ForwardIteratorB>
double manhattan_distance(
    ForwardIteratorA first_a, ForwardIteratorA last_a,
    ForwardIteratorB first_b, ForwardIteratorB last_b
) {
    return p_norm_distance( first_a, last_a, first_b, last_b, 1.0 );
}
 
inline double manhattan_distance(
    std::vector<double> const& vec_a, std::vector<double> const& vec_b
) {
    return p_norm_distance( vec_a.begin(), vec_a.end(), vec_b.begin(), vec_b.end(), 1.0 );
}
 
template <class ForwardIteratorA, class ForwardIteratorB>
double euclidean_distance(
    ForwardIteratorA first_a, ForwardIteratorA last_a,
    ForwardIteratorB first_b, ForwardIteratorB last_b
) {
    return p_norm_distance( first_a, last_a, first_b, last_b, 2.0 );
}
 
inline double euclidean_distance(
    std::vector<double> const& vec_a, std::vector<double> const& vec_b
) {
    return p_norm_distance( vec_a.begin(), vec_a.end(), vec_b.begin(), vec_b.end(), 2.0 );
}
 
template <class ForwardIteratorA, class ForwardIteratorB>
double maximum_distance(
    ForwardIteratorA first_a, ForwardIteratorA last_a,
    ForwardIteratorB first_b, ForwardIteratorB last_b
) {
    return p_norm_distance(
        first_a, last_a, first_b, last_b, std::numeric_limits<double>::infinity()
    );
}
 
inline double maximum_distance(
    std::vector<double> const& vec_a, std::vector<double> const& vec_b
) {
    return p_norm_distance(
        vec_a.begin(), vec_a.end(), vec_b.begin(), vec_b.end(),
        std::numeric_limits<double>::infinity()
    );
}
 
// =================================================================================================
//     Distances Matrices
// =================================================================================================
 
Matrix<double> p_norm_distance_matrix( Matrix<double> const& data, double p = 2.0 );
 
Matrix<double> manhattan_distance_matrix( Matrix<double> const& data );
 
Matrix<double> euclidean_distance_matrix( Matrix<double> const& data );
 
Matrix<double> maximum_distance_matrix( Matrix<double> const& data );
 
} // namespace utils
} // namespace genesis
 
#endif // include guard