logspace.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2006-2026 Knut Reinert & Freie Universität Berlin
// SPDX-FileCopyrightText: 2016-2026 Knut Reinert & MPI für molekulare Genetik
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <algorithm>
#include <cmath>
 
namespace raptor::logspace
{
 
constexpr double ln_2{0.693147180559945};
constexpr double negative_inf{-std::numeric_limits<double>::infinity()};
 
// Produces correct result if either term is -inf.
// Needs a check for the case where both terms are -inf.
[[nodiscard]] inline double add(double const log_x, double const log_y) noexcept
{
    double const max{std::max(log_x, log_y)};
    return max == negative_inf ? negative_inf : max + std::log1p(std::exp(-std::abs(log_x - log_y)));
}
 
template <typename... types>
[[nodiscard]] double add(double const log_x, double const log_y, types... logs) noexcept
{
    return add(add(log_x, log_y), logs...);
}
 
// log(x - y) = log_x + log(1 - std::exp(difference))
// expm1 is more accurate than using exp if the difference is close to 0.
// std::log1p(-std::exp(difference)) = std::log(1 + (-std::exp(difference)))
//                                   = std::log(1 - std::exp(difference))
// std::log(-std::expm1(difference)) = std::log(-(std::exp(difference) - 1))
//                                   = std::log(1 - std::exp(difference))
[[nodiscard]] inline double substract(double const log_x, double const log_y) noexcept
{
    double const difference{log_y - log_x};
    return log_x + (difference > -ln_2 ? std::log(-std::expm1(difference)) : std::log1p(-std::exp(difference)));
}
 
struct add_fn
{
    [[nodiscard]] double operator()(double const log_x, double const log_y) const noexcept
    {
        return add(log_x, log_y);
    };
};
 
} // namespace raptor::logspace