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-2018 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/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 ));

    double sum = 0.0;
    size_t cnt = 0;

    for_each_finite_pair( first_a, last_a, first_b, last_b, [&]( double val_a, double val_b ){
        if( std::isfinite( p )) {
            sum += std::pow( std::abs( val_a - val_b ), p );
        } else {
            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