#include <algorithm>
#include <boost/random.hpp>
#include <boost/bind.hpp>
#include <cstddef>
#include <functional>
#include <iostream>
#include <vector>


namespace my {


/**
 * \c Empirical discrete probability distribution.
 *
 * An empirical discrete distribution is an observed relative frequency
 * distribution that one intends to use to approximate the probabilities of a
 * random variable [1].
 * A discrete_distribution is modeled as a set of events \f$i\f$,
 * \f$0 \le i < n\f$, distributed according to the discrete probability function
 * \f[
 *   \Pr(i | p_0, \ldots , p_{n-1}) = p_i
 * \f]
 * Unless specified otherwise, the distribution parameters are calculated as:
 * \f[
 *   p_k = \frac{w_k}{S}, \quad for k = 0,\ldots, n-1
 * \f]
 * in which the values \f$w_k\f$ , commonly known as the \e weights, shall be
 * non-negative, non-NaN, and non-infinity.
 * Moreover, the following relation shall hold:
 * \f[
 *   0 < S = w_0 + \cdots + w_{n-1}
 * \f]
 *
 * References:
 * - [1] Carver et al. "Doing data analysis with SPSS version 16",
 *       Brooks/Cole, 4th edition, 2009
 *
 * \author Marco Guazzone, &lt;marco.guazzone@mfn.unipmn.it&gt;
 */
template <typename ResultT=int, typename WeightT=double>
class discrete_distribution
{

	public: typedef ResultT result_type;
	public: typedef WeightT weight_type;
	private: typedef ::std::vector<weight_type> probs_container;


	/**
 	 * \brief Constructs a discrete distribution object with \f$n = 1\f$ and
 	 * \f$p_0 = 1\f$.
 	 */
	public: discrete_distribution()
	{
		probs_.push_back(1);
		cum_probs_.push_back(1);
	}


	/// A constructor.
	public: template <typename WeightIteratorT>
		discrete_distribution(WeightIteratorT first_weight, WeightIteratorT last_weight)
	{
		weight_type weights_sum = 0;

		// Create PDF
		for (
			WeightIteratorT it = first_weight;
			it != last_weight;
			++it
		) {
			probs_.push_back(*it);
			weights_sum += *it;
		}

		// Normalize PDF and compute CDF
		weight_type cum_probs_sum = 0;
		for (
			typename probs_container::iterator it = probs_.begin();
			it != probs_.end();
			++it
		) {
			// Normalize probability
			*it /= weights_sum;

			// Compute CDF
			cum_probs_sum += *it;
			cum_probs_.push_back(cum_probs_sum);
		}
	}


	public: ::std::vector<weight_type> pdf() const
	{
		return probs_;
	}


	public: ::std::vector<weight_type> cdf() const
	{
		return cum_probs_;
	}


	// Compute the quantile value corresponding to the given probability value.
	public: result_type quantile(weight_type p) const
	{
		typename probs_container::const_iterator it;

		it = ::std::find_if(
			cum_probs_.begin(),
			cum_probs_.end(),
			::boost::bind(
				::std::greater_equal<weight_type>(),
				::_1,
				p
			)
		);

		if (it == cum_probs_.end())
		{
			return cum_probs_.size()-1;
		}

		return it-cum_probs_.begin();
	}


	/// Generate a random variate distributed according to this distribution and
	/// generated by means of the given uniform random number generator.
	public: template <typename UniformRandomGeneratorT>
		result_type operator()(UniformRandomGeneratorT& rng)
	{
#ifdef NDEBUG
		return quantile(::boost::uniform_01<weight_type>()(rng));
#else // NDEBUG
		weight_type r=::boost::uniform_01<weight_type>()(rng);
		::std::clog << "[my::discrete_distribution::rand] RNG: " << r << ::std::endl;
		return quantile(r);
#endif // NDEBUG
	}


	/// The probabilities container.
	private: probs_container probs_;
	/// The cumulative probabilities container.
	private: probs_container cum_probs_;

}; // discrete_distribution

} // Namespace my


int main()
{
	std::vector<double> freqs;
	freqs.push_back(1);
	freqs.push_back(4);
	freqs.push_back(5);

	boost::random::mt19937 rng;

	std::cout << "Testing Boost..." << std::endl;
	rng.seed(5489L);
	boost::random::discrete_distribution<> boost_dist(freqs.begin(), freqs.end());
	std::cout << "Random Event: " << boost_dist(rng) << std::endl;
	std::cout << "Random Event: " << boost_dist(rng) << std::endl;
	std::cout << "Random Event: " << boost_dist(rng) << std::endl;

	std::cout << std::endl;

	std::cout << "Testing My own version..." << std::endl;
	rng.seed(5489L);
	my::discrete_distribution<> my_dist(freqs.begin(), freqs.end());
	std::cout << "Random Event: " << my_dist(rng) << std::endl;
	std::cout << "Random Event: " << my_dist(rng) << std::endl;
	std::cout << "Random Event: " << my_dist(rng) << std::endl;

}