#pragma once

// #include <desc/distribution/FactorPruning.hxx>

#include <desc/distribution/IterativeScaling.hxx>
#include <desc/distribution/MaxEntFactor.hxx>
#include <desc/storage/Itemset.hxx>
// #include <relent/RelEnt.hxx>

#include <setcover/GreedySetcover.hxx>

#ifndef NDEBUG
#include <exception>
#endif

#if HAS_EXECUTION_POLICIES
#include <algorithm>
#include <execution>
#endif

namespace sd
{
namespace disc::mk2
{

template <typename Underlying_Factor_Type>
struct Factor
{
    using pattern_type = typename Underlying_Factor_Type::pattern_type;
    explicit Factor(size_t dim) : factor(dim) { range.reserve(dim); }
    itemset<pattern_type>  range;
    Underlying_Factor_Type factor;
};

template <typename U>
void join_factors(Factor<U>& f, const Factor<U>& g)
{
    for (auto& t : g.factor.itemsets.set) f.factor.insert_pattern(t.frequency, t.point);

    for (auto& t : g.factor.singletons.set) f.factor.insert_singleton(t.frequency, t.element);

    f.range.insert(g.range);
}

template <typename U>
struct Factorization
{
    using float_type = typename U::float_type;

    std::vector<Factor<U>> singleton_factors;
    std::vector<Factor<U>> factors;
};

template <typename U>
void init_singletons(std::vector<Factor<U>>& factors, size_t dim)
{
    factors.clear();
    factors.clear();
    factors.resize(dim, Factor<U>(dim));

    for (size_t i = 0; i < dim; ++i) {
        factors[i].range.insert(i);
        factors[i].factor.insert_singleton(0.5, i, true);
    }
}

template <typename U, typename float_type>
void insert_singleton(std::vector<Factor<U>>& factors,
                      size_t                  dim,
                      float_type              frequency,
                      const size_t            element,
                      bool                    estimate)
{
    bool add_new = true;
    for (auto& f : factors) {
        if (is_subset(element, f.range)) {
            f.factor.insert_singleton(frequency, element, estimate);
            if (is_singleton(f.range)) {
                add_new = false;
            }
        }
    }
    if (add_new) {
        auto& f = factors.emplace_back(dim);
        f.range.insert(element);
        f.factor.insert_singleton(frequency, element, estimate);
    }
}

template <typename U, typename F = double>
void estimate_model(std::vector<U>& phi, IterativeScalingSettings<F> const& opts = {})
{
#if HAS_EXECUTION_POLICIES
    std::for_each(std::execution::par_unseq, phi.begin(), phi.end(), [&](auto& f) {
        estimate_model(f.factor, opts);
    });
#else
#pragma omp parallel for
    for (size_t i = 0; i < factors.size(); ++i) {
        estimate_model(factors[i].factor, opts);
    }
#endif
}

template <typename U>
void erase_empty_factors(Factorization<U>& phi)
{
    phi.factors.erase(std::remove_if(phi.factors.begin(),
                                     phi.factors.end(),
                                     [](auto& f) { return f.range.empty(); }),
                      phi.factors.end());
    // phi.singleton_factors.erase(
    //     std::remove_if(phi.singleton_factors.begin(), phi.singleton_factors.end(), [](auto&
    //     f) { return f.range.empty(); }), phi.singleton_factors.end());
}

template <typename U, typename X, typename Visitor>
void factorize_classic(std::vector<Factor<U>>& phi, const X& x, Visitor&& f)
{
    for (size_t i = 0; i < phi.size(); ++i) {
        if (intersects(x, phi[i].range)) {
            f(phi[i], i);
        }
    }
}

template <typename U, typename X, typename Visitor>
void factorize_patterns_statically(const std::vector<Factor<U>>& factors, X& x, Visitor&& f)
{
    for (size_t i = 0; i < factors.size(); ++i) {
        if (intersects(x, factors[i].range)) {
            f(factors[i], i);
            setminus(x, factors[i].range);
        }
    }
}

template <typename Underlying_Factor_Type, typename X, typename Visitor>
void factorize_statically(const Factorization<Underlying_Factor_Type>& phi,
                          const X&                                     x,
                          Visitor&&                                    f)
{
    using pattern_type = typename Underlying_Factor_Type::pattern_type;

    thread_local itemset<pattern_type> y;
    y.clear();
    y.assign(x);

    factorize_patterns_statically(
        phi.factors, y, [&](const auto& phi_i, size_t i) { f(phi_i, i, false); });
    for_each_item(y, [&](size_t i) { f(phi.singleton_factors[i], i, true); });
}

template <typename U, typename V, typename Underlying_Factor_Type = MaxEntFactor<U, V>>
struct StaticFactorModel
{
    using float_type            = V;
    using pattern_type          = U;
    using index_type            = std::size_t;
    using underlying_model_type = Underlying_Factor_Type;
    using factor_type           = Factor<underlying_model_type>;

    constexpr static bool enable_factor_pruning = false;

    explicit StaticFactorModel(size_t dim) : dim(dim) { init(dim); }
    StaticFactorModel() = default;

    size_t dim              = 0;
    size_t max_factor_size  = 5;
    size_t max_factor_width = 8;

    Factorization<Underlying_Factor_Type> phi;

    size_t num_itemsets() const { return phi.factors.size(); }
    size_t dimension() const { return dim; }
    size_t size() const { return num_itemsets() + dimension(); }

    void init(size_t d)
    {
        dim = d;
        clear();
        init_singletons(phi.singleton_factors, dim);
    }

    void clear()
    {
        phi.factors.clear();
        phi.singleton_factors.clear();
    }

    void insert_singleton(float_type frequency, const index_type element, bool estimate)
    {
        disc::mk2::insert_singleton(phi.singleton_factors, dim, frequency, element, estimate);
    }

    template <typename T>
    void insert_singleton(float_type frequency, const T& t, bool estimate)
    {
        insert_singleton(frequency, static_cast<index_type>(front(t)), estimate);
    }

    void prune_factor_if(factor_type& f, size_t max_factor_size)
    {
        if constexpr (enable_factor_pruning) {
            if (f.factor.itemsets.set.size() > 3) {
                prune_factor(f, max_factor_size);
            }
        }
    }

    template <typename T>
    void insert_pattern(float_type frequency,
                        const T&   t,
                        size_t     max_factor_size,
                        size_t     max_factor_width,
                        bool       estimate)
    {

        std::vector<std::pair<size_t, bool>> selection;
        selection.reserve(count(t));

        bool found_superset = false;
        factorize_statically(phi, t, [&](const auto& f, size_t i, bool from_singleton) {
            if (is_subset(t, f.range) && !from_singleton) {
                if (f.factor.itemsets.set.size() < max_factor_size) {
                    auto& phi_i = phi.singleton_factors[i];
                    if (estimate) prune_factor_if(phi_i, max_factor_size);
                    phi_i.factor.insert(frequency, t, estimate);
                }
                found_superset = true;
                return;
            }

            if (!found_superset) selection.emplace_back(i, from_singleton);
        });

        if (selection.empty() || found_superset) {
            return;
        }

        factor_type next(dim);

        for (const auto & [ i, s ] : selection) {
            join_factors(next, s ? phi.singleton_factors[i] : phi.factors[i]);
        }

        if (count(next.range) > max_factor_width ||
            next.factor.itemsets.set.size() > max_factor_size)
        {
            // #ifndef NDEBUG
            //             throw std::domain_error{"pattern too large or factor is full"};
            // #endif
            return;
        }

        if (estimate) prune_factor_if(next, max_factor_size);
        next.factor.insert(frequency, t, estimate);

        for (auto[i, s] : selection) {
            if (!s) phi.factors[i].range.clear();
        }
        selection.clear();

        erase_empty_factors(phi);

        phi.factors.emplace_back(std::move(next));
    }

    template <typename T>
    void insert_pattern(float_type frequency, const T& t, bool estimate)
    {
        insert_pattern(frequency, t, max_factor_size, max_factor_width, estimate);
    }

    template <typename T>
    void insert(float_type frequency, const T& t, bool estimate)
    {
        if (is_singleton(t))
            insert_singleton(frequency, t, estimate);
        else
            insert_pattern(frequency, t, estimate);
    }

    template <typename T>
    bool is_allowed(const T& t, size_t max_factor_size, size_t max_factor_width) const
    {
        size_t total_size     = 0;
        size_t total_width    = 0;
        bool   found_superset = false;

        factorize_statically(phi, t, [&](const auto& f, size_t, bool) {
            if (is_subset(t, f.range)) {
                found_superset = true;
                total_size     = f.factor.itemsets.set.size();
                total_width    = f.factor.singletons.set.size();
            }
            else if (!found_superset)
            {
                total_size += f.factor.itemsets.set.size();
                total_width += count(f.range);
            }
        });

        return total_size < max_factor_size && total_width < max_factor_width;
    }

    template <typename T>
    bool is_allowed(const T& t) const
    {
        return is_allowed(t, max_factor_size, max_factor_width);
    }
};

template <typename S, typename T, typename U, typename X, typename Visitor>
void factorize(const StaticFactorModel<S, T, U>& pr, const X& x, Visitor&& f)
{
    factorize_statically(pr.phi, x, std::forward<Visitor>(f));
}
template <typename U, typename V, typename X, typename Visitor>
void factorize_patterns(const StaticFactorModel<U, V>& pr, const X& x, Visitor&& f)
{
    thread_local itemset<typename StaticFactorModel<U, V>::pattern_type> y;
    y.clear();
    y.assign(x);

    factorize_patterns_statically(pr.phi.factors, y, f);
}
template <typename S, typename T, typename U, typename F = double>
void estimate_model(StaticFactorModel<S, T, U>& m, IterativeScalingSettings<F> const& opts = {})
{
    estimate_model(m.phi.factors, opts);
}

template <typename U, typename V, typename Underlying_Factor_Type = MaxEntFactor<U, V>>
struct ClassicStaticFactorModel
{
    using float_type            = V;
    using value_type            = V;
    using pattern_type          = U;
    using index_type            = std::size_t;
    using underlying_model_type = Underlying_Factor_Type;
    // using cache_type            = Cache<U, V>;
    using factor_type = Factor<underlying_model_type>;

    constexpr static bool enable_factor_pruning = false;

    explicit ClassicStaticFactorModel(size_t dim) : dim(dim) { init(dim); }
    ClassicStaticFactorModel() = default;

    size_t dim              = 0;
    size_t count_sets       = 0;
    size_t max_factor_size  = 5;
    size_t max_factor_width = 8;

    std::vector<factor_type> factors;

    size_t num_itemsets() const { return count_sets; }
    size_t dimension() const { return dim; }
    size_t size() const { return num_itemsets() + dimension(); }

    void init(size_t dim)
    {
        this->dim = dim;
        clear();
        init_singletons(factors, dim);
    }

    void clear() { factors.clear(); }

    void insert_singleton(value_type frequency, const index_type element, bool estimate)
    {
        disc::mk2::insert_singleton(factors, dim, frequency, element, estimate);
    }

    template <typename T>
    void insert_singleton(value_type frequency, const T& t, bool estimate)
    {
        insert_singleton(frequency, static_cast<index_type>(front(t)), estimate);
    }

    void prune_factor_if(factor_type& f, size_t max_factor_size)
    {
        if constexpr (enable_factor_pruning) {
            if (f.factor.itemsets.set.size() > 3) {
                prune_factor(f, max_factor_size);
            }
        }
    }

    template <typename T>
    void insert_pattern(value_type frequency,
                        const T&   t,
                        size_t     max_factor_size,
                        size_t     max_factor_width,
                        bool       estimate)
    {
        std::vector<size_t> selection;
        selection.reserve(count(t));

        for (size_t i = 0, n = factors.size(); i < n; ++i) {
            auto& f = factors[i];
            if (is_subset(t, f.range)) {
                if (f.factor.itemsets.set.size() < max_factor_size) {
                    if (estimate) prune_factor_if(f, max_factor_size);
                    f.factor.insert(frequency, t, estimate);
                }
                return;
            }
            if (intersects(t, f.range)) selection.push_back(i);
        }

        if (selection.empty()) {
            return;
        }

        factor_type next(dim);

        for (const auto& s : selection) {
            join_factors(next, factors[s]);
        }

        if (count(next.range) > max_factor_width ||
            next.factor.itemsets.set.size() > max_factor_size)
        {
            // #ifndef NDEBUG
            //             throw std::domain_error{"pattern too large or factor is full"};
            // #endif
            return;
        }

        if (estimate) prune_factor_if(next, max_factor_size);
        next.factor.insert(frequency, t, estimate);

        for (auto i : selection) {
            factors[i].range.clear();
        }
        selection.clear();

        erase_empty_factors();

        factors.emplace_back(std::move(next));
        ++count_sets;
    }

    template <typename T>
    void insert_pattern(value_type frequency, const T& t, bool estimate)
    {
        insert_pattern(frequency, t, max_factor_size, max_factor_width, estimate);
    }

    template <typename T>
    void insert(value_type frequency, const T& t, bool estimate)
    {
        if (is_singleton(t))
            insert_singleton(frequency, t, estimate);
        else
            insert_pattern(frequency, t, estimate);
    }

    template <typename T>
    bool is_allowed(const T& t, size_t max_factor_size, size_t max_factor_width) const
    {
        size_t total_size  = 0;
        size_t total_width = 0;
        for (auto& f : factors) {
            if (is_subset(t, f.range)) {
                return f.factor.itemsets.set.size() < max_factor_size;
            }
            if (intersects(t, f.range)) {
                // this works, because all factors are covering disjoint sets
                total_size += f.factor.itemsets.set.size();
                total_width += count(f.range);
                if (total_size >= max_factor_size || total_width >= max_factor_width)
                    return false;
            }
        }

        return true;
    }

    template <typename T>
    bool is_allowed(const T& t) const
    {
        return is_allowed(t, max_factor_size, max_factor_width);
    }

    underlying_model_type as_single_factor() const
    {
        factor_type next(dim);
        for (const auto& s : factors) {
            join_factors(next, s);
        }
        return next.factor;
    }

    void erase_empty_factors()
    {
        factors.erase(std::remove_if(factors.begin(),
                                     factors.end(),
                                     [](auto& f) { return f.range.empty(); }),
                      factors.end());
        // for (size_t i = 0; i < factors.size();)
        // {
        //     if (factors[i].range.empty())
        //     {
        //         if (i != factors.size() - 1)
        //         {
        //             std::swap(factors[i], factors.back());
        //         }

        //         factors.pop_back();
        //     }
        //     else
        //     {
        //         ++i;
        //     }
        // }
    }
};
template <typename S, typename T, typename U, typename X, typename Visitor>
void factorize(const ClassicStaticFactorModel<S, T, U>& pr, const X& x, Visitor&& v)
{
    factorize_classically(pr.factors, x, [&](const auto& f, size_t i) {
        v(f, i, f.factor.singletons.set.size() == 1);
    });
}
template <typename U, typename V, typename X, typename Visitor>
void factorize_patterns(const ClassicStaticFactorModel<U, V>& pr, const X& x, Visitor&& v)
{
    factorize_classically(pr.factors, x, [&](const auto& f, size_t i) {
        if (f.factor.singletons.set.size() != 1) v(f, i);
    });
}
template <typename S, typename T, typename U, typename F = double>
void estimate_model(ClassicStaticFactorModel<S, T, U>& m,
                    IterativeScalingSettings<F> const& opts = {})
{
    estimate_model(m.factors, opts);
}

struct SetcoverWeightFunction
{
    template <typename Set, typename To_Cover>
    float operator()(size_t, const Set& f, const To_Cover& y) const noexcept
    {
        const auto s = size_of_intersection(y, f.range);

        if (s > 1 && f.factor.itemsets.get_precomputed_expectation(y)) {
            return std::numeric_limits<float>::max();
        }
        const auto t = f.factor.itemsets.width();
        return s == t ? (static_cast<float>(s)) : -1.f;
    }
};

struct SuperSetcoverWeightFunction
{
    template <typename Set, typename To_Cover>
    float operator()(size_t, const Set& f, const To_Cover& y) const noexcept
    {
        const auto s = size_of_intersection(y, f.range);

        if (s > 1 && f.factor.itemsets.get_precomputed_expectation(y)) {
            return std::numeric_limits<float>::max();
        }

        return static_cast<float>(s);
    }
};

struct WeightedSuperSetcoverWeightFunction
{
    template <typename Set, typename To_Cover>
    float operator()(size_t, const Set& f, const To_Cover& y) const noexcept
    {
        const auto s = size_of_intersection(y, f.range);

        if (s > 1 && f.factor.itemsets.get_precomputed_expectation(y)) {
            return std::numeric_limits<float>::max();
        }

        const auto t = f.factor.itemsets.width();
        return static_cast<float>(s) / (t - s + 1);
    }
};

struct SetMinusFactor
{
    template <typename To_Cover, typename Set>
    void operator()(To_Cover& s, const Set& f) const noexcept
    {
        sd::setminus(s, f.range);
    }
};

template <typename U, typename X, typename Visitor>
void factorize_patterns_dynamically(const std::vector<Factor<U>>& factors, X& x, Visitor&& f)
{
    using W = SuperSetcoverWeightFunction;
    setcover::greedy_set_cover(factors,
                               x,
                               [&](size_t i) { f(factors[i], i, false); },
                               W{},
                               SetMinusFactor{});
}

template <typename U, typename X, typename Visitor>
void factorize_dynamically(const Factorization<U>& phi, const X& x, Visitor&& f)
{
    thread_local itemset<typename U::pattern_type> y;
    y.clear();
    y.assign(x);

    factorize_patterns_dynamically(phi.factors, y, f);
    for_each_item(y, [&](size_t i) { f(phi.singleton_factors[i], i, true); });
}

template <typename factor_type>
void erase_uncovered_singletons(factor_type& f, bool estimate)
{
    auto&       s   = f.factor.singletons.set;
    const auto& r   = f.range;
    const auto  pos = std::remove_if(
        s.begin(), s.end(), [&r](const auto& i) { return !is_subset(i.element, r); });

    if (pos != s.end()) {
        s.erase(pos, s.end());

        if (estimate) estimate_model(f.factor);
        f.factor.itemsets.num_singletons = f.factor.singletons.size();
    }
}

template <typename factor_type, typename float_type, typename T>
void create_mindiv_factor(factor_type& next,
                          float_type   frequency,
                          const T&     t,
                          size_t       max_factor_size,
                          size_t       max_factor_width,
                          bool         estimate)
{
    assert(!is_singleton(t));
    // singletons + most informative next pattern

    factor_type replacement(next.factor.singletons.dim);

    replacement.factor.singletons              = next.factor.singletons;
    replacement.factor.itemsets.num_singletons = replacement.factor.singletons.size();
    replacement.range.insert(t);
    replacement.factor.insert(frequency, t, estimate);

    auto& xs = next.factor.itemsets.set;

    thread_local itemset<tag_dense> in_use;
    in_use.clear();

    // greedy algorithm to select the sets from next that happen to increase the likelihood
    for (size_t count = 1, n = xs.size(); count <= max_factor_size;) {
        std::pair best{n, float_type(0)};
        for (size_t j = 0; j < n; ++j) {
            if (in_use.test(j)) {
                continue;
            }

            if (!intersects(xs[j].point, t)) {
                in_use.insert(j);
                continue;
            }

            if (size_of_union(xs[j].point, replacement.range) > max_factor_width) {
                in_use.insert(j);
                continue;
            }
            using std::log2;

            auto p = expectation(replacement.factor, xs[j].point);
            auto q = xs[j].frequency;
            auto g = (q * log2(q / p)) + (p * log2(p / q)); // assignment_score

            if (g >= best.second) {
                best = {j, g};
            }
        }

        if (best.first < n && best.second != 0) {
            ++count;
            in_use.insert(best.first);
            const auto& y = xs[best.first].point;
            const auto& q = xs[best.first].frequency;
            replacement.range.insert(y);
            replacement.factor.insert(q, y, estimate);
        }
        else
        {
            break;
        }
    }

    mk2::erase_uncovered_singletons(replacement, estimate);
    next = std::move(replacement);
}

template <typename factor_type, typename T>
void prune_factor_set_cover_based(factor_type& next,
                                  const T&     t,
                                  size_t       max_factor_size,
                                  bool         estimate)
{
    auto&       is = next.factor.itemsets.set;
    factor_type replacement(next.factor.singletons.dim);

    auto weight = [](size_t, const auto& xs, const auto& ys) {
        size_t s = sd::size_of_intersection(xs.point, ys);
        return s;
        // return s == count(xs.point) ? s : size_t(0);
    };
    auto setminus = [](auto& ys, const auto& xs) { sd::setminus(ys, xs.point); };
    auto report   = [&](size_t i) {
        replacement.factor.insert(is[i].frequency, is[i].point);
        replacement.range.insert(is[i].point);
    };

    // setcover::greedy_set_cover_by_value(is, t, report, {max_factor_size - 1}, weight,
    // setminus);
    disc::itemset<disc::tag_dense> ignore(is.size());

    setcover::greedy_set_cover_no_minus(ignore, is, t, report, {max_factor_size - 1}, weight);

    for (const auto& i : next.factor.singletons.set) {
        if (is_subset(i.element, replacement.range)) {
            replacement.factor.insert_singleton(i.frequency, i.element);
        }
    }
    if (estimate) estimate_model(replacement.factor);
    next = std::move(replacement);
}

template <typename factor_type, typename T>
void prune_factor_overlapping_based(factor_type& next,
                                    const T&     t,
                                    size_t       max_factor_size,
                                    bool         estimate)
{
    // keep most essential sets
    auto& is = next.factor.itemsets.set;
    std::sort(std::begin(is), std::end(is), [&t](const auto& a, const auto& b) {
        return size_of_intersection(a.point, t) > size_of_intersection(b.point, t);
    });
    is.resize(max_factor_size);
    next.range.clear();
    for (const auto& i : is) {
        next.range.insert(i.point);
    }
    next.range.insert(t);
    auto& sset = next.factor.singletons.set;
    sset.erase(std::remove_if(std::begin(sset),
                              std::end(sset),
                              [&](const auto& s) { return !is_subset(s.element, next.range); }),
               std::end(sset));
    next.factor.itemsets.num_singletons = sset.size();

    if (estimate) estimate_model(next.factor);
}

template <typename U, typename V, typename Underlying_Factor_Type = MaxEntFactor<U, V>>
struct RelaxedFactorModel
{
    using float_type            = V;
    using pattern_type          = U;
    using index_type            = std::size_t;
    using underlying_model_type = Underlying_Factor_Type;
    using factor_type           = Factor<underlying_model_type>;

    explicit RelaxedFactorModel(size_t dim) : dim(dim) { init(dim); }
    RelaxedFactorModel() = default;

    size_t dim              = 0;
    size_t max_factor_size  = 5;
    size_t max_factor_width = 8;

    Factorization<Underlying_Factor_Type> phi;

    size_t num_itemsets() const { return phi.factors.size(); }
    size_t dimension() const { return dim; }
    size_t size() const { return num_itemsets() + dimension(); }

    void init(size_t d)
    {
        dim = d;
        clear();
        init_singletons(phi.singleton_factors, dim);
    }

    void clear()
    {
        phi.factors.clear();
        phi.singleton_factors.clear();
    }

    void insert_singleton(float_type frequency, const index_type element, bool estimate)
    {
        disc::mk2::insert_singleton(phi.singleton_factors, dim, frequency, element, estimate);
    }

    template <typename T>
    void insert_singleton(float_type frequency, const T& t, bool estimate)
    {
        insert_singleton(frequency, static_cast<index_type>(front(t)), estimate);
    }

    template <typename T>
    void insert_pattern_join_only_singletons(float_type frequency, const T& t, bool estimate)
    {
        factor_type next(this->dim);
        for_each_item(t, [&](size_t i) { join_factors(next, phi.singleton_factors[i]); });

        next.range.insert(t);
        next.factor.insert(frequency, t, estimate);
        phi.factors.emplace_back(std::move(next));
    }

    template <typename T>
    void insert_pattern(float_type fr,
                        const T&   t,
                        size_t     max_factor_size,
                        size_t     max_factor_width,
                        bool       estimate)
    {

        if (max_factor_size <= 1) {
            insert_pattern_join_only_singletons(fr, t, estimate);
            return;
        }

        // itemset<U>  x_buf(t);
        factor_type next(this->dim);

        // factorize_dynamically(
        //     *this, this->dim, x_buf, [&](auto& f, size_t, bool) { join_factors(next, f); });

        factorize_dynamically(
            phi, t, [&](const auto& f, size_t, bool) { join_factors(next, f); });

        const auto next_elems = next.factor.itemsets.set.size();
        const auto next_items = next.factor.singletons.set.size();

        if ((next_elems >= max_factor_size || next_items > max_factor_width)) {
            mk2::create_mindiv_factor(next, fr, t, max_factor_size, max_factor_width, true);
        }
        else
        {
            next.factor.insert(fr, t, estimate);
            next.range.insert(t);
        }

        this->phi.factors.emplace_back(std::move(next));
    }

    template <typename T>
    void insert_pattern(float_type frequency, const T& t, bool estimate)
    {
        insert_pattern(frequency, t, max_factor_size, max_factor_width, estimate);
    }

    template <typename T>
    void insert(float_type frequency, const T& t, bool estimate)
    {
        if (is_singleton(t))
            insert_singleton(frequency, t, estimate);
        else
            insert_pattern(frequency, t, estimate);
    }

    template <typename T>
    bool is_allowed(const T&, size_t, size_t) const
    {
        return true;
    }

    template <typename T>
    bool is_allowed(const T& t) const
    {
        return is_allowed(t, max_factor_size, max_factor_width);
    }
};
template <typename S, typename T, typename U, typename X, typename Visitor>
void factorize(const RelaxedFactorModel<S, T, U>& pr, const X& x, Visitor&& f)
{
    factorize_dynamically(pr.phi, x, std::forward<Visitor>(f));
}
template <typename U, typename V, typename X, typename Visitor>
void factorize_patterns(const RelaxedFactorModel<U, V>& pr, const X& x, Visitor&& f)
{
    thread_local itemset<typename RelaxedFactorModel<U, V>::pattern_type> y;
    y.clear();
    y.assign(x);

    factorize_patterns_dynamically(pr.phi.factors, y, f);
}
template <typename S, typename T, typename U, typename F = double>
void estimate_model(RelaxedFactorModel<S, T, U>&       m,
                    IterativeScalingSettings<F> const& opts = {})
{
    estimate_model(m.phi.factors, opts);
}

template <typename Model, typename X>
auto log_expectation(Model const& m, X const& x)
{
    using std::log2;

    using float_type   = typename Model::float_type;
    using pattern_type = typename Model::pattern_type;

    thread_local itemset<pattern_type> part;

    float_type fr = 0;
    factorize(m, x, [&](const auto& f, size_t, bool s) {
        if (s) {
            fr += log2(f.factor.singletons.set.front().probability);
        }
        else
        {
            part.clear();
            intersection(x, f.range, part);
            fr += log2(expectation(f.factor, part));
        }
    });
    return fr;
}

template <typename Model, typename X>
auto log_probability(Model const& m, X const& x)
{
    using std::log2;

    using float_type   = typename Model::float_type;
    using pattern_type = typename Model::pattern_type;

    thread_local itemset<pattern_type> part;

    float_type fr = 0;
    factorize(m, x, [&](const auto& f, size_t, bool s) {
        if (s) {
            fr += log2(f.factor.singletons.set.front().probability);
        }
        else
        {
            part.clear();
            intersection(x, f.range, part);
            fr += log2(probability(f.factor, part));
        }
    });
    return fr;
}

template <typename U, typename X, typename Visitor>
void factorize_singletons(const Factorization<U>& phi, const X& x, Visitor&& v)
{
    for_each_item(x, [&](size_t i) { v(phi.singleton_factors[i], i); });
}

template <typename S, typename T, typename U, typename X, typename Visitor>
void factorize_singletons(const RelaxedFactorModel<S, T, U>& pr, const X& x, Visitor&& v)
{
    factorize_singletons(pr.phi, x, std::forward<Visitor>(v));
}
template <typename S, typename T, typename U, typename X, typename Visitor>
void factorize_singletons(const StaticFactorModel<S, T, U>& pr, const X& x, Visitor&& v)
{
    factorize_singletons(pr.phi, x, std::forward<Visitor>(v));
}

template <typename U, typename X>
auto log_probability_of_absent_items(const std::vector<Factor<U>>& phi, const X& x)
{
    using std::log2;
    using float_type = typename U::float_type;
    float_type p     = 0;
    for (size_t i = 0; i < phi.size(); ++i) {
        const auto& e = phi[i].factor.singletons.set;
        if (e.size() == 1 && !is_subset(e.front().element, x)) {
            p += log2(1 - e.front().probability);
        }
    }
    return p;
}

template <typename U, typename X>
auto log_probability_of_absent_items(const Factorization<U>& phi, const X& x)
{
    using std::log2;
    using float_type = typename U::float_type;
    float_type p     = 0;
    for (size_t i = 0; i < phi.singleton_factors.size(); ++i) {
        const auto& e = phi.singleton_factors[i].factor.singletons.set;
        if (!is_subset(e.front().element, x)) {
            p += log2(1 - e.front().probability);
        }
    }
    return p;
}

template <typename S, typename T, typename U, typename X>
auto log_probability_of_absent_items(const RelaxedFactorModel<S, T, U>& pr, const X& x)
{
    return log_probability_of_absent_items(pr.phi, x);
}
template <typename S, typename T, typename U, typename X>
auto log_probability_of_absent_items(const StaticFactorModel<S, T, U>& pr, const X& x)
{
    return log_probability_of_absent_items(pr.phi, x);
}
template <typename S, typename T, typename U, typename X>
auto log_probability_of_absent_items(const ClassicStaticFactorModel<S, T, U>& pr, const X& x)
{
    return log_probability_of_absent_items(pr.factors, x);
}

template <typename Model, typename X>
auto log_expectation_generalized_set(const Model& m, const X& x)
{
    using std::log2;
    return log_expectation(m, x) + log_probability_of_absent_items(m, x);
}

template <typename Model, typename X>
auto expectation(Model const& m, X const& x)
{
    using std::exp2;
    return exp2(log_expectation(m, x));
}

template <typename Model, typename X>
auto probability(Model const& m, X const& x)
{
    using std::exp2;
    return exp2(log_probability(m, x));
}
template <typename Model, typename X>
auto expectation_generalized_set(Model const& m, X const& x)
{
    using std::exp2;
    return exp2(log_expectation_generalized_set(m, x));
}

} // namespace disc::mk2
} // namespace sd