Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // extended_p_square.hpp
Chris@16: //
Chris@16: //  Copyright 2005 Daniel Egloff. Distributed under the Boost
Chris@16: //  Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_HPP_DE_01_01_2006
Chris@16: #define BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_HPP_DE_01_01_2006
Chris@16: 
Chris@16: #include <vector>
Chris@16: #include <functional>
Chris@16: #include <boost/range/begin.hpp>
Chris@16: #include <boost/range/end.hpp>
Chris@16: #include <boost/range/iterator_range.hpp>
Chris@16: #include <boost/iterator/transform_iterator.hpp>
Chris@16: #include <boost/iterator/counting_iterator.hpp>
Chris@16: #include <boost/iterator/permutation_iterator.hpp>
Chris@16: #include <boost/parameter/keyword.hpp>
Chris@16: #include <boost/mpl/placeholders.hpp>
Chris@16: #include <boost/accumulators/accumulators_fwd.hpp>
Chris@16: #include <boost/accumulators/framework/extractor.hpp>
Chris@16: #include <boost/accumulators/numeric/functional.hpp>
Chris@16: #include <boost/accumulators/framework/parameters/sample.hpp>
Chris@16: #include <boost/accumulators/framework/depends_on.hpp>
Chris@16: #include <boost/accumulators/statistics_fwd.hpp>
Chris@16: #include <boost/accumulators/statistics/count.hpp>
Chris@16: #include <boost/accumulators/statistics/times2_iterator.hpp>
Chris@16: 
Chris@16: namespace boost { namespace accumulators
Chris@16: {
Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // probabilities named parameter
Chris@16: //
Chris@16: BOOST_PARAMETER_NESTED_KEYWORD(tag, extended_p_square_probabilities, probabilities)
Chris@16: 
Chris@16: BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_probabilities)
Chris@16: 
Chris@16: namespace impl
Chris@16: {
Chris@16:     ///////////////////////////////////////////////////////////////////////////////
Chris@16:     // extended_p_square_impl
Chris@16:     //  multiple quantile estimation
Chris@16:     /**
Chris@16:         @brief Multiple quantile estimation with the extended \f$P^2\f$ algorithm
Chris@16: 
Chris@16:         Extended \f$P^2\f$ algorithm for estimation of several quantiles without storing samples.
Chris@16:         Assume that \f$m\f$ quantiles \f$\xi_{p_1}, \ldots, \xi_{p_m}\f$ are to be estimated.
Chris@16:         Instead of storing the whole sample cumulative distribution, the algorithm maintains only
Chris@16:         \f$m+2\f$ principal markers and \f$m+1\f$ middle markers, whose positions are updated
Chris@16:         with each sample and whose heights are adjusted (if necessary) using a piecewise-parablic
Chris@16:         formula. The heights of these central markers are the current estimates of the quantiles
Chris@16:         and returned as an iterator range.
Chris@16: 
Chris@16:         For further details, see
Chris@16: 
Chris@16:         K. E. E. Raatikainen, Simultaneous estimation of several quantiles, Simulation, Volume 49,
Chris@16:         Number 4 (October), 1986, p. 159-164.
Chris@16: 
Chris@16:         The extended \f$ P^2 \f$ algorithm generalizes the \f$ P^2 \f$ algorithm of
Chris@16: 
Chris@16:         R. Jain and I. Chlamtac, The P^2 algorithm for dynamic calculation of quantiles and
Chris@16:         histograms without storing observations, Communications of the ACM,
Chris@16:         Volume 28 (October), Number 10, 1985, p. 1076-1085.
Chris@16: 
Chris@16:         @param extended_p_square_probabilities A vector of quantile probabilities.
Chris@16:     */
Chris@16:     template<typename Sample>
Chris@16:     struct extended_p_square_impl
Chris@16:       : accumulator_base
Chris@16:     {
Chris@16:         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
Chris@16:         typedef std::vector<float_type> array_type;
Chris@16:         // for boost::result_of
Chris@16:         typedef iterator_range<
Chris@16:             detail::lvalue_index_iterator<
Chris@16:                 permutation_iterator<
Chris@16:                     typename array_type::const_iterator
Chris@16:                   , detail::times2_iterator
Chris@16:                 >
Chris@16:             >
Chris@16:         > result_type;
Chris@16: 
Chris@16:         template<typename Args>
Chris@16:         extended_p_square_impl(Args const &args)
Chris@16:           : probabilities(
Chris@16:                 boost::begin(args[extended_p_square_probabilities])
Chris@16:               , boost::end(args[extended_p_square_probabilities])
Chris@16:             )
Chris@16:           , heights(2 * probabilities.size() + 3)
Chris@16:           , actual_positions(heights.size())
Chris@16:           , desired_positions(heights.size())
Chris@16:           , positions_increments(heights.size())
Chris@16:         {
Chris@16:             std::size_t num_quantiles = this->probabilities.size();
Chris@16:             std::size_t num_markers = this->heights.size();
Chris@16: 
Chris@16:             for(std::size_t i = 0; i < num_markers; ++i)
Chris@16:             {
Chris@16:                 this->actual_positions[i] = i + 1;
Chris@16:             }
Chris@16: 
Chris@16:             this->positions_increments[0] = 0.;
Chris@16:             this->positions_increments[num_markers - 1] = 1.;
Chris@16: 
Chris@16:             for(std::size_t i = 0; i < num_quantiles; ++i)
Chris@16:             {
Chris@16:                 this->positions_increments[2 * i + 2] = probabilities[i];
Chris@16:             }
Chris@16: 
Chris@16:             for(std::size_t i = 0; i <= num_quantiles; ++i)
Chris@16:             {
Chris@16:                 this->positions_increments[2 * i + 1] =
Chris@16:                     0.5 * (this->positions_increments[2 * i] + this->positions_increments[2 * i + 2]);
Chris@16:             }
Chris@16: 
Chris@16:             for(std::size_t i = 0; i < num_markers; ++i)
Chris@16:             {
Chris@16:                 this->desired_positions[i] = 1. + 2. * (num_quantiles + 1.) * this->positions_increments[i];
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         template<typename Args>
Chris@16:         void operator ()(Args const &args)
Chris@16:         {
Chris@16:             std::size_t cnt = count(args);
Chris@16: 
Chris@16:             // m+2 principal markers and m+1 middle markers
Chris@16:             std::size_t num_markers = 2 * this->probabilities.size() + 3;
Chris@16: 
Chris@16:             // first accumulate num_markers samples
Chris@16:             if(cnt <= num_markers)
Chris@16:             {
Chris@16:                 this->heights[cnt - 1] = args[sample];
Chris@16: 
Chris@16:                 // complete the initialization of heights by sorting
Chris@16:                 if(cnt == num_markers)
Chris@16:                 {
Chris@16:                     std::sort(this->heights.begin(), this->heights.end());
Chris@16:                 }
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                 std::size_t sample_cell = 1;
Chris@16: 
Chris@16:                 // find cell k = sample_cell such that heights[k-1] <= sample < heights[k]
Chris@16:                 if(args[sample] < this->heights[0])
Chris@16:                 {
Chris@16:                     this->heights[0] = args[sample];
Chris@16:                     sample_cell = 1;
Chris@16:                 }
Chris@16:                 else if(args[sample] >= this->heights[num_markers - 1])
Chris@16:                 {
Chris@16:                     this->heights[num_markers - 1] = args[sample];
Chris@16:                     sample_cell = num_markers - 1;
Chris@16:                 }
Chris@16:                 else
Chris@16:                 {
Chris@16:                     typedef typename array_type::iterator iterator;
Chris@16:                     iterator it = std::upper_bound(
Chris@16:                         this->heights.begin()
Chris@16:                       , this->heights.end()
Chris@16:                       , args[sample]
Chris@16:                     );
Chris@16: 
Chris@16:                     sample_cell = std::distance(this->heights.begin(), it);
Chris@16:                 }
Chris@16: 
Chris@16:                 // update actual positions of all markers above sample_cell index
Chris@16:                 for(std::size_t i = sample_cell; i < num_markers; ++i)
Chris@16:                 {
Chris@16:                     ++this->actual_positions[i];
Chris@16:                 }
Chris@16: 
Chris@16:                 // update desired positions of all markers
Chris@16:                 for(std::size_t i = 0; i < num_markers; ++i)
Chris@16:                 {
Chris@16:                     this->desired_positions[i] += this->positions_increments[i];
Chris@16:                 }
Chris@16: 
Chris@16:                 // adjust heights and actual positions of markers 1 to num_markers-2 if necessary
Chris@16:                 for(std::size_t i = 1; i <= num_markers - 2; ++i)
Chris@16:                 {
Chris@16:                     // offset to desired position
Chris@16:                     float_type d = this->desired_positions[i] - this->actual_positions[i];
Chris@16: 
Chris@16:                     // offset to next position
Chris@16:                     float_type dp = this->actual_positions[i+1] - this->actual_positions[i];
Chris@16: 
Chris@16:                     // offset to previous position
Chris@16:                     float_type dm = this->actual_positions[i-1] - this->actual_positions[i];
Chris@16: 
Chris@16:                     // height ds
Chris@16:                     float_type hp = (this->heights[i+1] - this->heights[i]) / dp;
Chris@16:                     float_type hm = (this->heights[i-1] - this->heights[i]) / dm;
Chris@16: 
Chris@16:                     if((d >= 1 && dp > 1) || (d <= -1 && dm < -1))
Chris@16:                     {
Chris@16:                         short sign_d = static_cast<short>(d / std::abs(d));
Chris@16: 
Chris@16:                         float_type h = this->heights[i] + sign_d / (dp - dm) * ((sign_d - dm)*hp
Chris@16:                                      + (dp - sign_d) * hm);
Chris@16: 
Chris@16:                         // try adjusting heights[i] using p-squared formula
Chris@16:                         if(this->heights[i - 1] < h && h < this->heights[i + 1])
Chris@16:                         {
Chris@16:                             this->heights[i] = h;
Chris@16:                         }
Chris@16:                         else
Chris@16:                         {
Chris@16:                             // use linear formula
Chris@16:                             if(d > 0)
Chris@16:                             {
Chris@16:                                 this->heights[i] += hp;
Chris@16:                             }
Chris@16:                             if(d < 0)
Chris@16:                             {
Chris@16:                                 this->heights[i] -= hm;
Chris@16:                             }
Chris@16:                         }
Chris@16:                         this->actual_positions[i] += sign_d;
Chris@16:                     }
Chris@16:                 }
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         result_type result(dont_care) const
Chris@16:         {
Chris@16:             // for i in [1,probabilities.size()], return heights[i * 2]
Chris@16:             detail::times2_iterator idx_begin = detail::make_times2_iterator(1);
Chris@16:             detail::times2_iterator idx_end = detail::make_times2_iterator(this->probabilities.size() + 1);
Chris@16: 
Chris@16:             return result_type(
Chris@16:                 make_permutation_iterator(this->heights.begin(), idx_begin)
Chris@16:               , make_permutation_iterator(this->heights.begin(), idx_end)
Chris@16:             );
Chris@16:         }
Chris@16: 
Chris@16:     private:
Chris@16:         array_type probabilities;         // the quantile probabilities
Chris@16:         array_type heights;               // q_i
Chris@16:         array_type actual_positions;      // n_i
Chris@16:         array_type desired_positions;     // d_i
Chris@16:         array_type positions_increments;  // f_i
Chris@16:     };
Chris@16: 
Chris@16: } // namespace impl
Chris@16: 
Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // tag::extended_p_square
Chris@16: //
Chris@16: namespace tag
Chris@16: {
Chris@16:     struct extended_p_square
Chris@16:       : depends_on<count>
Chris@16:       , extended_p_square_probabilities
Chris@16:     {
Chris@16:         typedef accumulators::impl::extended_p_square_impl<mpl::_1> impl;
Chris@16: 
Chris@16:         #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED
Chris@16:         /// tag::extended_p_square::probabilities named parameter
Chris@16:         static boost::parameter::keyword<tag::probabilities> const probabilities;
Chris@16:         #endif
Chris@16:     };
Chris@16: }
Chris@16: 
Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // extract::extended_p_square
Chris@16: //
Chris@16: namespace extract
Chris@16: {
Chris@16:     extractor<tag::extended_p_square> const extended_p_square = {};
Chris@16: 
Chris@16:     BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square)
Chris@16: }
Chris@16: 
Chris@16: using extract::extended_p_square;
Chris@16: 
Chris@16: // So that extended_p_square can be automatically substituted with
Chris@16: // weighted_extended_p_square when the weight parameter is non-void
Chris@16: template<>
Chris@16: struct as_weighted_feature<tag::extended_p_square>
Chris@16: {
Chris@16:     typedef tag::weighted_extended_p_square type;
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct feature_of<tag::weighted_extended_p_square>
Chris@16:   : feature_of<tag::extended_p_square>
Chris@16: {
Chris@16: };
Chris@16: 
Chris@16: }} // namespace boost::accumulators
Chris@16: 
Chris@16: #endif