///////////////////////////////////////////////////////////////////////////////
// joint_histogram.hpp
//
//  Copyright 2007 Hirotaka Niitsuma. 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_JOINT_HISTOGRAM_HPP_DE_01_01_2007
#define BOOST_ACCUMULATORS_STATISTICS_JOINT_HISTOGRAM_HPP_DE_01_01_2007

#include <vector>
#include <deque>
#include <limits>
#include <numeric>
#include <functional>
#include <complex>
#include <boost/mpl/assert.hpp>
#include <boost/mpl/bool.hpp>
#include <boost/range.hpp>
#include <boost/parameter/keyword.hpp>
#include <boost/mpl/placeholders.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/type_traits/is_scalar.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/accumulators/framework/accumulator_base.hpp>
#include <boost/accumulators/framework/extractor.hpp>
#include <boost/numeric/functional.hpp>
#include <boost/accumulators/framework/parameters/sample.hpp>
#include <boost/accumulators/statistics_fwd.hpp>
#include <boost/accumulators/statistics/count.hpp>
#if USE_ACC_MIN_MAX_IN_ACC_DENSITY
#include <boost/accumulators/statistics/max.hpp>
#include <boost/accumulators/statistics/min.hpp>
#endif
#include <boost/accumulators/statistics/sample_cache.hpp>
#include <boost/accumulators/statistics/density.hpp>

#define FOR_TEST_NIITSUMA 0
namespace boost { namespace accumulators
{

namespace impl
{
    ///////////////////////////////////////////////////////////////////////////////
    // joint_histogram_impl
    //
     template<typename Sample, typename VariateType, typename VariateTag>
    struct joint_histogram_impl
      : accumulator_base
    {

        typedef typename numeric::functional::average<Sample, std::size_t>::result_type sample_type;
		typedef std::deque<sample_type> cache_type;
		typedef iterator_range<typename cache_type::iterator> range_cache_type;

		typedef std::vector<std::pair<sample_type, sample_type> > histogram_type;
		typedef iterator_range<typename histogram_type::iterator> range_histogram_type;

		typedef boost::numeric::ublas::matrix<sample_type>  result_type;
			

  
        template<typename Args>
        joint_histogram_impl(Args const &args)
          : 
		
		cache_size(args[sample_cache_size])
		,num_bins(args[density_num_bins])
		,joint_histogram( boost::numeric::ublas::zero_matrix<sample_type> ( num_bins + 2,num_bins + 2))

        {}
        

		int value_to_index_of_histogram_range(range_histogram_type &hist_range, sample_type &sample_data )
		{				
				   
                   if (sample_data < hist_range[1].first)	
                    {     						
						return 0;
                    }
                    else if (sample_data >=  hist_range[this->num_bins + 1].first)
                    {                       
					   return num_bins + 1;
                    }    
                    else
                    {
						for (typename range_histogram_type::const_iterator hit=hist_range.begin() 
							; hit != hist_range.end(); ++hit )
						{
							if( sample_data  < hit->first )
							{
								return ( hit-hist_range.begin()-1 ) ;
							}
								
						}
                    }			
		}

        template<typename Args>
        void operator ()(Args const &args)
        {

			std::size_t cnt = count(args);
		
            if (cnt == this->cache_size)
            {

				this->cache_of_sample		= sample_cache(args);
				this->density_of_sample		= density(args);											 

				extractor<tag::sample_cache_of_variates<VariateType, VariateTag> > const some_sample_cache_of_variates = {};
				extractor<tag::density_of_variates<VariateType, VariateTag> > const some_density_of_variates = {};				
				this->cache_of_variates		= some_sample_cache_of_variates(args);
				this->density_of_variates	= some_density_of_variates(args);

				range_cache_type::const_iterator iter				= this->cache_of_sample.begin();
				range_cache_type::const_iterator iter_of_variates	= this->cache_of_variates.begin();
				int index=0,index_of_variates=0;

                for (; iter != this->cache_of_sample.end(); ++iter,++iter_of_variates )
                {
					index				= value_to_index_of_histogram_range(this->density_of_sample		, *iter );
					index_of_variates	= value_to_index_of_histogram_range(this->density_of_variates	, *iter_of_variates);
					joint_histogram(index,  index_of_variates )++;
                }
            }
            // Add each subsequent sample to the correct bin
            else if (cnt > this->cache_size)
            {

				sample_type sample_data=args[sample];
				sample_type sample_variates = args[parameter::keyword<VariateTag>::get()];
				int index=0,index_of_variates=0;
					index				= value_to_index_of_histogram_range(this->density_of_sample		, sample_data );
					index_of_variates	= value_to_index_of_histogram_range(this->density_of_variates	, sample_variates);

					joint_histogram(index,  index_of_variates )++;
            }       
        }

        result_type result(dont_care) const
        {
			return this->joint_histogram;
        }

    private:
		std::size_t				num_bins;        // number of bins
		std::size_t				cache_size;      // number of cached samples		
	
		range_cache_type		cache_of_sample;
		range_cache_type		cache_of_variates;
		range_histogram_type	density_of_sample;
		range_histogram_type	density_of_variates;
		result_type				joint_histogram;		
		
    };
    
} // namespace impl

///////////////////////////////////////////////////////////////////////////////
// tag::joint_histogram
//
namespace tag
{
    template<typename VariateType, typename VariateTag>
    struct joint_histogram
      : depends_on<
			count 
#if USE_ACC_MIN_MAX_IN_ACC_DENSITY
			,min,max
			,min_of_variates<VariateType, VariateTag> 
			,max_of_variates<VariateType, VariateTag>
#endif
			,sample_cache
			,density			
			,sample_cache_of_variates<VariateType, VariateTag>
			,density_of_variates<VariateType, VariateTag>		   
	  >
    {
        typedef accumulators::impl::joint_histogram_impl<mpl::_1, VariateType, VariateTag> impl;
    };  
    
    struct abstract_joint_histogram 
      : depends_on<>
    {
    };
}

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

using extract::joint_histogram;

template<typename VariateType, typename VariateTag>
struct feature_of<tag::joint_histogram<VariateType, VariateTag> >
  : feature_of<tag::abstract_joint_histogram>
{
};

}} // namespace boost::accumulators

#endif