distance_metrics.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>
#include <stdexcept>
#include <iterator>
#include <vector>
 
namespace statcpp {
 
template <typename Iterator1, typename Iterator2>
double euclidean_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::euclidean_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::euclidean_distance: empty sequences");
    }
 
    double sum_sq = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double diff = static_cast<double>(*it1) - static_cast<double>(*it2);
        sum_sq += diff * diff;
        ++it1;
        ++it2;
    }
 
    return std::sqrt(sum_sq);
}
 
template <typename Iterator1, typename Iterator2, typename Proj1, typename Proj2>
double euclidean_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2,
                          Proj1 proj1, Proj2 proj2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::euclidean_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::euclidean_distance: empty sequences");
    }
 
    double sum_sq = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double v1 = static_cast<double>(proj1(*it1));
        double v2 = static_cast<double>(proj2(*it2));
        double diff = v1 - v2;
        sum_sq += diff * diff;
        ++it1;
        ++it2;
    }
 
    return std::sqrt(sum_sq);
}
 
template <typename Iterator1, typename Iterator2>
double manhattan_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::manhattan_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::manhattan_distance: empty sequences");
    }
 
    double sum_abs = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double diff = static_cast<double>(*it1) - static_cast<double>(*it2);
        sum_abs += std::abs(diff);
        ++it1;
        ++it2;
    }
 
    return sum_abs;
}
 
template <typename Iterator1, typename Iterator2, typename Proj1, typename Proj2>
double manhattan_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2,
                          Proj1 proj1, Proj2 proj2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::manhattan_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::manhattan_distance: empty sequences");
    }
 
    double sum_abs = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double v1 = static_cast<double>(proj1(*it1));
        double v2 = static_cast<double>(proj2(*it2));
        double diff = v1 - v2;
        sum_abs += std::abs(diff);
        ++it1;
        ++it2;
    }
 
    return sum_abs;
}
 
template <typename Iterator1, typename Iterator2>
double cosine_similarity(Iterator1 first1, Iterator1 last1,
                         Iterator2 first2, Iterator2 last2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::cosine_similarity: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::cosine_similarity: empty sequences");
    }
 
    double dot_product = 0.0;
    double norm1_sq = 0.0;
    double norm2_sq = 0.0;
 
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double v1 = static_cast<double>(*it1);
        double v2 = static_cast<double>(*it2);
 
        dot_product += v1 * v2;
        norm1_sq += v1 * v1;
        norm2_sq += v2 * v2;
 
        ++it1;
        ++it2;
    }
 
    if (norm1_sq == 0.0 || norm2_sq == 0.0) {
        throw std::invalid_argument("statcpp::cosine_similarity: zero vector encountered");
    }
 
    return dot_product / (std::sqrt(norm1_sq) * std::sqrt(norm2_sq));
}
 
template <typename Iterator1, typename Iterator2, typename Proj1, typename Proj2>
double cosine_similarity(Iterator1 first1, Iterator1 last1,
                         Iterator2 first2, Iterator2 last2,
                         Proj1 proj1, Proj2 proj2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::cosine_similarity: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::cosine_similarity: empty sequences");
    }
 
    double dot_product = 0.0;
    double norm1_sq = 0.0;
    double norm2_sq = 0.0;
 
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double v1 = static_cast<double>(proj1(*it1));
        double v2 = static_cast<double>(proj2(*it2));
 
        dot_product += v1 * v2;
        norm1_sq += v1 * v1;
        norm2_sq += v2 * v2;
 
        ++it1;
        ++it2;
    }
 
    if (norm1_sq == 0.0 || norm2_sq == 0.0) {
        throw std::invalid_argument("statcpp::cosine_similarity: zero vector encountered");
    }
 
    return dot_product / (std::sqrt(norm1_sq) * std::sqrt(norm2_sq));
}
 
template <typename Iterator1, typename Iterator2>
double cosine_distance(Iterator1 first1, Iterator1 last1,
                       Iterator2 first2, Iterator2 last2)
{
    return 1.0 - cosine_similarity(first1, last1, first2, last2);
}
 
template <typename Iterator1, typename Iterator2, typename Proj1, typename Proj2>
double cosine_distance(Iterator1 first1, Iterator1 last1,
                       Iterator2 first2, Iterator2 last2,
                       Proj1 proj1, Proj2 proj2)
{
    return 1.0 - cosine_similarity(first1, last1, first2, last2, proj1, proj2);
}
 
inline double mahalanobis_distance(const std::vector<double>& x,
                                   const std::vector<double>& mean,
                                   const std::vector<std::vector<double>>& cov_matrix)
{
    if (x.size() != mean.size()) {
        throw std::invalid_argument("statcpp::mahalanobis_distance: x and mean must have same dimension");
    }
    if (x.size() != 2) {
        throw std::invalid_argument("statcpp::mahalanobis_distance: only 2D is currently supported");
    }
    if (cov_matrix.size() != 2 || cov_matrix[0].size() != 2 || cov_matrix[1].size() != 2) {
        throw std::invalid_argument("statcpp::mahalanobis_distance: covariance matrix must be 2x2");
    }
 
    // Compute difference vector: d = x - mean
    double d0 = x[0] - mean[0];
    double d1 = x[1] - mean[1];
 
    // Compute inverse of 2x2 covariance matrix
    // Σ = [a b]
    //     [c d]
    // Σ⁻¹ = (1/det) * [d  -b]
    //                 [-c  a]
 
    double a = cov_matrix[0][0];
    double b = cov_matrix[0][1];
    double c = cov_matrix[1][0];
    double d = cov_matrix[1][1];
 
    double det = a * d - b * c;
 
    if (std::abs(det) < 1e-10) {
        throw std::invalid_argument("statcpp::mahalanobis_distance: covariance matrix is singular");
    }
 
    // Inverse matrix elements
    double inv_a = d / det;
    double inv_b = -b / det;
    double inv_c = -c / det;
    double inv_d = a / det;
 
    // Compute dᵀ Σ⁻¹ d
    // First compute Σ⁻¹ d
    double temp0 = inv_a * d0 + inv_b * d1;
    double temp1 = inv_c * d0 + inv_d * d1;
 
    // Then compute dᵀ (Σ⁻¹ d)
    double mahalanobis_sq = d0 * temp0 + d1 * temp1;
 
    if (mahalanobis_sq < 0.0) {
        // Due to numerical errors, ensure non-negative
        mahalanobis_sq = 0.0;
    }
 
    return std::sqrt(mahalanobis_sq);
}
 
template <typename Iterator1, typename Iterator2>
double minkowski_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2,
                          double p)
{
    if (p < 1.0) {
        throw std::invalid_argument("statcpp::minkowski_distance: p must be >= 1");
    }
 
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::minkowski_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::minkowski_distance: empty sequences");
    }
 
    double sum = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double diff = std::abs(static_cast<double>(*it1) - static_cast<double>(*it2));
        sum += std::pow(diff, p);
        ++it1;
        ++it2;
    }
 
    return std::pow(sum, 1.0 / p);
}
 
template <typename Iterator1, typename Iterator2, typename Proj1, typename Proj2>
double minkowski_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2,
                          double p, Proj1 proj1, Proj2 proj2)
{
    if (p < 1.0) {
        throw std::invalid_argument("statcpp::minkowski_distance: p must be >= 1");
    }
 
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::minkowski_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::minkowski_distance: empty sequences");
    }
 
    double sum = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double v1 = static_cast<double>(proj1(*it1));
        double v2 = static_cast<double>(proj2(*it2));
        double diff = std::abs(v1 - v2);
        sum += std::pow(diff, p);
        ++it1;
        ++it2;
    }
 
    return std::pow(sum, 1.0 / p);
}
 
template <typename Iterator1, typename Iterator2>
double chebyshev_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::chebyshev_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::chebyshev_distance: empty sequences");
    }
 
    double max_diff = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double diff = std::abs(static_cast<double>(*it1) - static_cast<double>(*it2));
        if (diff > max_diff) {
            max_diff = diff;
        }
        ++it1;
        ++it2;
    }
 
    return max_diff;
}
 
template <typename Iterator1, typename Iterator2, typename Proj1, typename Proj2>
double chebyshev_distance(Iterator1 first1, Iterator1 last1,
                          Iterator2 first2, Iterator2 last2,
                          Proj1 proj1, Proj2 proj2)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        throw std::invalid_argument("statcpp::chebyshev_distance: sequences must have the same length");
    }
    if (n1 == 0) {
        throw std::invalid_argument("statcpp::chebyshev_distance: empty sequences");
    }
 
    double max_diff = 0.0;
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double v1 = static_cast<double>(proj1(*it1));
        double v2 = static_cast<double>(proj2(*it2));
        double diff = std::abs(v1 - v2);
        if (diff > max_diff) {
            max_diff = diff;
        }
        ++it1;
        ++it2;
    }
 
    return max_diff;
}
 
} // namespace statcpp