
#pragma once

#include <disc/disc/Composition.hxx>
#include <disc/disc/PatternsetResult.hxx>
#include <disc/distribution/Distribution.hxx>
#include <disc/storage/Dataset.hxx>
#include <disc/utilities/Support.hxx>
#include <disc/utilities/UniversalIntEncoding.hxx>

#include <math/nchoosek.hxx>

#include <algorithm>
#include <cmath>
#include <cstddef>

namespace sd::disc::bic
{

template <typename Trait>
auto encode_model_bic(const Composition<Trait>& c) //_aaai20
{
    // n is constant, add or remove if from score to your desire.
    // times 2 because we consider theta and assignment of patterns as dof

    const auto n = c.data.size();
    const auto s = c.summary.size();
    const auto k = c.data.num_components();
    const auto d = c.data.dim;
    const auto m = s - d;

    using float_type    = typename Trait::float_type;
    const float_type df = k * (2 * m + d) + n;
    // const float_type df = k * m * 2; // + d + n const

    return std::log2(n) * df / 2.;
}

template <typename Trait>
auto encode_model_bic_new(const Composition<Trait>& c)
{
    const auto n  = c.data.size();
    const auto k  = c.data.num_components();
    const auto d  = c.data.dim;
    const auto s  = c.summary.size();
    const auto m  = s - d;
    const auto tr = trace(c.assignment) - d * k; // singletons are free

    using float_type    = typename Trait::float_type;
    const float_type df = tr + m * k; // + n

    return std::log2(n) * df / 2.;
}

template <typename Trait>
auto encode_model_bic(const PatternsetResult<Trait>& c)
{
    return std::log2(c.data.size()) * (c.summary.size() - c.data.dim) / 2.;
}

#if 0

template <typename C>
auto additional_cost_bic(const C& c ) -> typename C::float_type
{
    using data_type = std::decay_t<decltype(c.data)>;

    size_t df = 1;
    if constexpr (meta::has_components_member_fn<data_type>::value)
    {
        auto k = c.data.num_components();
        df += k > 1 ? k : 0;
    }
    return std::log2(c.data.size()) * df / 2;
}
#endif

} // namespace sd::disc::bic