ProbabilityDistribution.hpp

// Copyright 2021 Arnaud Becheler    <abechele@umich.edu>
 
/************************************************************************ * This program is free software; you can
 *redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free
 *Software Foundation; either version 2 of the License, or    * (at your option) any later version. *
 *                                                                      *
 ***************************************************************************/
 
#ifndef __OCCUPANCY_SPECTRUM_DISTRIBUTION_H_INCLUDED__
#define __OCCUPANCY_SPECTRUM_DISTRIBUTION_H_INCLUDED__
 
#include "Support.hpp"
#include "editor_policy.hpp"
#include "filter_policy.hpp"
 
#include <assert.h>
#include <cmath>    // std::pow
#include <iostream> // << operator implementation
#include <random>   // discrete_distribution
#include <utility>  // std::forward, std::move
#include <vector>
 
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/random.hpp> // boost::multiprecision::pow
 
namespace quetzal
{
namespace coalescence
{
namespace occupancy_spectrum
{
template <class FilterPolicy = filter_policy::return_always_true, class EditorPolicy = editor_policy::identity>
class ProbabilityDistribution
{
  private:
    using BigInt = boost::multiprecision::cpp_int;
    using BigFloat = boost::multiprecision::cpp_dec_float_50;
    using self_type = ProbabilityDistribution<FilterPolicy, EditorPolicy>;
    using spectrum_type = OccupancySpectrum;
    using support_type = std::vector<spectrum_type>;
    using probabilities_type = std::vector<double>;
    int m_k;
    int m_N;
    FilterPolicy m_pred;
    support_type m_support;
    probabilities_type m_probas;
    mutable std::discrete_distribution<size_t> m_distribution;
 
  public:
    ProbabilityDistribution() = default;
    ProbabilityDistribution(int n, int m, FilterPolicy pred = FilterPolicy()) : m_k(n), m_N(m), m_pred(pred)
    {
        auto callback = make_callback(m_k, m_N);
        Support s(m_k, m_N);
        s.generate(callback);
        m_distribution = std::discrete_distribution<size_t>(m_probas.cbegin(), m_probas.cend());
        assert(m_support.size() == m_probas.size());
    }
    ProbabilityDistribution(const self_type &) = delete;
    ProbabilityDistribution(self_type &&) = default;
    self_type &operator=(self_type &&) = default;
    self_type &operator=(const self_type &) = delete;
    template <typename Generator> const spectrum_type &operator()(Generator &g) const
    {
        return m_support[m_distribution(g)];
    }
    int n() const
    {
        return m_k;
    }
    int m() const
    {
        return m_N;
    }
    const support_type &support() const
    {
        return m_support;
    }
    const probabilities_type &weights() const
    {
        return m_probas;
    }
    static double compute_probability(int n, int m, spectrum_type const &m_js)
    {
        BigFloat numerator(Factorial(m) * Factorial(n));
        BigFloat denominator(1.0);
        int j = 0;
        for (auto m_j : m_js)
        {
            denominator *= boost::multiprecision::pow(BigFloat(Factorial(j)), BigFloat(m_j)) * BigFloat(Factorial(m_j));
            ++j;
        }
        denominator *= BigFloat(std::pow(m, n));
        BigFloat result = numerator / denominator;
        return static_cast<double>(result);
    }
    friend std::ostream &operator<<(std::ostream &os, self_type const &dist)
    {
        for (int i = 0; i < dist.support().size(); ++i)
        {
            // os << "P( \t" << dist.support()[i] << "\t) = " << dist.weights()[i] << "\n";
            auto v = dist.support()[i];
            os << "P( \t";
            if (!v.empty())
            {
                os << '[';
                std::copy(v.begin(), v.end(), std::ostream_iterator<int>(os, ", "));
                os << "\b\b]";
            }
            os << "\t) = " << dist.weights()[i] << "\n";
        }
        return os;
    }
 
  private:
    auto make_callback(int k, int N)
    {
        return [this, k, N](spectrum_type &&M_j) {
            auto p = compute_probability(k, N, M_j);
            if (this->m_pred(p))
            {
                this->m_support.push_back(EditorPolicy::edit(std::move(M_j)));
                this->m_probas.push_back(p);
            }
        };
    }
    static BigInt Factorial(int num)
    {
        BigInt fact = 1;
        for (int i = num; i > 1; --i)
        {
            fact *= i;
        }
        return fact;
    }
};
} // namespace occupancy_spectrum
} // namespace coalescence
} // namespace quetzal
 
#endif