///////////////////////////////////////////////////////////////////////////////
// simple_density.hpp
//
//  Copyright 2006 Daniel Egloff, Olivier Gygi. Distributed under the Boost
//  Modified - 2008 Milosz Marian Hulboj. Distributed under the Boost
//  Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_ACCUMULATORS_STATISTICS_SIMPLE_DENSITY_HPP_CH_27_01_2008
#define BOOST_ACCUMULATORS_STATISTICS_SIMPLE_DENSITY_HPP_CH_27_01_2008

#include <vector>
#include <limits>
#include <functional>
#include <boost/range.hpp>
#include <boost/parameter/keyword.hpp>
#include <boost/mpl/placeholders.hpp>
#include <boost/accumulators/framework/accumulator_base.hpp>
#include <boost/accumulators/framework/extractor.hpp>
#include <boost/accumulators/numeric/functional.hpp>
#include <boost/accumulators/framework/parameters/sample.hpp>
#include <boost/accumulators/framework/depends_on.hpp>
#include <boost/accumulators/statistics_fwd.hpp>
#include <boost/accumulators/statistics/count.hpp>

namespace boost { namespace accumulators
{

///////////////////////////////////////////////////////////////////////////////
// cache_size and num_bins named parameters
//
BOOST_PARAMETER_NESTED_KEYWORD(tag, simple_density_lower_bound, lower_bound)
BOOST_PARAMETER_NESTED_KEYWORD(tag, simple_density_upper_bound, upper_bound)
BOOST_PARAMETER_NESTED_KEYWORD(tag, simple_density_num_bins, num_bins)

namespace impl
{
    ///////////////////////////////////////////////////////////////////////////////
    // simple_density_impl
    //  simple_density histogram
    /**
        @brief Histogram density estimator

        The histogram density estimator returns a histogram of the sample distribution. The positions and sizes of the bins
        are predefined by the user. The range between the minimum and the maximum provided by user  is subdivided into a 
        specifiable number of bins (num_bins) of same size. Additionally, an under- and an overflow bin is added to capture
        future under- and overflow samples. 
        At the end, a range of std::pair is being returned, where each pair contains the position of the bin (lower bound) 
        and the samples count (normalized with the total number of samples).
        Note:
        $x\in B_{i}\Leftrightarrow\left\{ \begin{array}{c}x\geq B_{i}^{min}\\x<B_{i}^{max}\end{array}\right.$
        From which follows that the samples equal to provided upper bound will fall into the overflow bin.

        @param  simple_density_lower_bound Lower bound of the first bin
        @param  simple_density_upper_bound Upper bound of the last bin
        @param  simple_density_num_bins Number of bins (two additional bins collect under- and overflow samples).
    */
    template<typename Sample>
    struct simple_density_impl
      : accumulator_base
    {
        typedef typename numeric::functional::average<Sample, std::size_t>::result_type float_type;
        typedef std::vector<std::pair<float_type, float_type> > histogram_type;
        typedef std::vector<float_type> array_type;
        // for boost::result_of
        typedef iterator_range<typename histogram_type::iterator> result_type;

        template<typename Args>
        simple_density_impl(Args const &args)
            : num_bins(args[simple_density_num_bins])
            , minimum(args[simple_density_lower_bound])
            , maximum(args[simple_density_upper_bound])
            , samples_in_bin(num_bins + 2, 0.)
            , bin_positions(num_bins + 2)
            , histogram(
                num_bins + 2
              , std::make_pair(
                    numeric::average(args[sample | Sample()],(std::size_t)1)
                  , numeric::average(args[sample | Sample()],(std::size_t)1)
                )
              )
            , is_dirty(true)
        {
            // Create num_bins bins of the same size between the minimum and maximum and two
            // additional bins for over- and under- flows. Prepare the lower bounds.
            float_type bin_size = numeric::average(maximum - minimum, this->num_bins );
            
            // determine bin positions (their lower bounds)
            for (std::size_t i = 0; i < this->num_bins + 2; ++i)
            {
                this->bin_positions[i] = minimum + (i - 1.) * bin_size;
            }
        }

        template<typename Args>
        void operator ()(Args const &args)
        {
            this->is_dirty = true;

            std::size_t cnt = count(args);

            if (args[sample] < this->bin_positions[1])
            {
                ++(this->samples_in_bin[0]);
            }
            else if (args[sample] >= this->bin_positions[this->num_bins + 1])
            {
                ++(this->samples_in_bin[this->num_bins + 1]);
            }
            else
            {
                typename array_type::iterator it = std::upper_bound(
                    this->bin_positions.begin()
                    , this->bin_positions.end()
                    , args[sample]
                );

                std::size_t d = std::distance(this->bin_positions.begin(), it);
                ++(this->samples_in_bin[d - 1]);
            }
        }

        template<typename Args>
        result_type result(Args const &args) const
        {
            if (this->is_dirty)
            {
                this->is_dirty = false;

                // creates a vector of std::pair where each pair i holds
                // the values bin_positions[i] (x-axis of histogram) and
                // samples_in_bin[i] / cnt (y-axis of histogram).

                for (std::size_t i = 0; i < this->num_bins + 2; ++i)
                {
                    this->histogram[i] = std::make_pair(this->bin_positions[i], numeric::average(this->samples_in_bin[i], count(args)));
                }
            }
            // returns a range of pairs
            return make_iterator_range(this->histogram);
        }

    private:
        std::size_t            num_bins;        // number of bins
        float_type             minimum;         // Lower bound of first bin
        float_type             maximum;         // Upper bound of last bin == lower bound for overflow bin
        array_type             samples_in_bin;  // number of samples in each bin
        array_type             bin_positions;   // lower bounds of bins
        mutable histogram_type histogram;       // histogram
        mutable bool is_dirty;
    };

} // namespace impl

///////////////////////////////////////////////////////////////////////////////
// tag::simple_density
//
namespace tag
{
    struct simple_density
      : depends_on<count>
      , simple_density_lower_bound
      , simple_density_upper_bound
      , simple_density_num_bins
    {
        /// INTERNAL ONLY
        ///
        typedef accumulators::impl::simple_density_impl<mpl::_1> impl;

        #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED
        /// tag::density::lower_bound named parameter
        /// tag::density::upper_bound named parameter
        /// tag::density::num_bins named parameter
        static boost::parameter::keyword<simple_density_lower_bound> const lower_bound;
        static boost::parameter::keyword<simple_density_upper_bound> const upper_bound;
        static boost::parameter::keyword<simple_density_num_bins> const num_bins;
        #endif
    };
}

///////////////////////////////////////////////////////////////////////////////
// extract::simple_density
//
namespace extract
{
    extractor<tag::simple_density> const simple_density = {};
}

using extract::simple_density;

}} // namespace boost::accumulators

#endif