#pragma once
#include <elisa/Binomial.hxx>
#include <math/nchoosek.hxx>

#include <cmath>

namespace sd::disc::fwer
{

size_t isum_binomial_coeff(size_t N, size_t k)
{
    size_t nck = 1, sum = 1;
    for (size_t i = 1; i <= k; ++i) {
        nck *= (static_cast<double>(N) - i + 1) / i;
        sum += nck;
    }
    return sum;
}

template <typename T>
auto log_alpha_bonferroni_correction(size_t, size_t N)
{
    return - T(N);
}
template <typename T>
auto log_alpha_webb_correction(size_t k, size_t N)
{
    return - T(k) - log_nchoosek(N, k);
}
template <typename T>
auto log_alpha_isumbc_correction(size_t K, size_t N)
{
    return - log(isum_binomial_coeff(N, K));
}
template <typename T>
auto log_alpha_geo_correction(size_t K, size_t N)
{
    using std::log;
    const T k = static_cast<T>(K);

    /*
        const T n = static_cast<T>(N);
        auto lugo = log(alpha) -log_nchoosek(n, k) + log(n - (k - 1)) - log(n - (2 * k - 1));
        auto halevi = log(alpha) -log_nchoosek(N, k + 1) + log(k + 1) - log(N - 2 * k - 1);
    */

    // log(a / (e^k (N/k)^k))
    return - k * (log(N) - log(k) + 1);
}

template <typename T>
auto log_alpha_correction(const T& alpha, size_t K, size_t N)
{
    using std::log;
    // return log_alpha_bonferroni_correction(alpha, K, N);
    // return log_alpha_webb_correction(alpha, K, N);
    // return log_alpha_isumbc_correction(alpha, K, N);
    return log(alpha) + log_alpha_geo_correction<T>(K, N);
}

} // namespace sd::disc