accumulators/statistics/p_square_cumul_dist.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // p_square_cumulative_distribution.hpp
0003 //
0004 //  Copyright 2005 Daniel Egloff, Olivier Gygi. Distributed under the Boost
0005 //  Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_CUMUL_DIST_HPP_DE_01_01_2006
0009 #define BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_CUMUL_DIST_HPP_DE_01_01_2006
0010
0011 #include <vector>
0012 #include <functional>
0013 #include <boost/parameter/keyword.hpp>
0014 #include <boost/range.hpp>
0015 #include <boost/mpl/placeholders.hpp>
0016 #include <boost/accumulators/accumulators_fwd.hpp>
0017 #include <boost/accumulators/framework/accumulator_base.hpp>
0018 #include <boost/accumulators/framework/extractor.hpp>
0019 #include <boost/accumulators/numeric/functional.hpp>
0020 #include <boost/accumulators/framework/parameters/sample.hpp>
0021 #include <boost/accumulators/statistics_fwd.hpp>
0022 #include <boost/accumulators/statistics/count.hpp>
0023 #include <boost/serialization/vector.hpp>
0024 #include <boost/serialization/utility.hpp>
0025
0026 namespace boost { namespace accumulators
0027 {
0028 ///////////////////////////////////////////////////////////////////////////////
0029 // num_cells named parameter
0030 //
0031 BOOST_PARAMETER_NESTED_KEYWORD(tag, p_square_cumulative_distribution_num_cells, num_cells)
0032
0033 BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_cumulative_distribution_num_cells)
0034
0035 namespace impl
0036 {
0037     ///////////////////////////////////////////////////////////////////////////////
0038     // p_square_cumulative_distribution_impl
0039     //  cumulative_distribution calculation (as histogram)
0040     /**
0041         @brief Histogram calculation of the cumulative distribution with the \f$P^2\f$ algorithm
0042
0043         A histogram of the sample cumulative distribution is computed dynamically without storing samples
0044         based on the \f$ P^2 \f$ algorithm. The returned histogram has a specifiable amount (num_cells)
0045         equiprobable (and not equal-sized) cells.
0046
0047         For further details, see
0048
0049         R. Jain and I. Chlamtac, The P^2 algorithm for dynamic calculation of quantiles and
0050         histograms without storing observations, Communications of the ACM,
0051         Volume 28 (October), Number 10, 1985, p. 1076-1085.
0052
0053         @param p_square_cumulative_distribution_num_cells.
0054     */
0055     template<typename Sample>
0056     struct p_square_cumulative_distribution_impl
0057       : accumulator_base
0058     {
0059         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
0060         typedef std::vector<float_type> array_type;
0061         typedef std::vector<std::pair<float_type, float_type> > histogram_type;
0062         // for boost::result_of
0063         typedef iterator_range<typename histogram_type::iterator> result_type;
0064
0065         template<typename Args>
0066         p_square_cumulative_distribution_impl(Args const &args)
0067           : num_cells(args[p_square_cumulative_distribution_num_cells])
0068           , heights(num_cells + 1)
0069           , actual_positions(num_cells + 1)
0070           , desired_positions(num_cells + 1)
0071           , positions_increments(num_cells + 1)
0072           , histogram(num_cells + 1)
0073           , is_dirty(true)
0074         {
0075             std::size_t b = this->num_cells;
0076
0077             for (std::size_t i = 0; i < b + 1; ++i)
0078             {
0079                 this->actual_positions[i] = i + 1.;
0080                 this->desired_positions[i] = i + 1.;
0081                 this->positions_increments[i] = numeric::fdiv(i, b);
0082             }
0083         }
0084
0085         template<typename Args>
0086         void operator ()(Args const &args)
0087         {
0088             this->is_dirty = true;
0089
0090             std::size_t cnt = count(args);
0091             std::size_t sample_cell = 1; // k
0092             std::size_t b = this->num_cells;
0093
0094             // accumulate num_cells + 1 first samples
0095             if (cnt <= b + 1)
0096             {
0097                 this->heights[cnt - 1] = args[sample];
0098
0099                 // complete the initialization of heights by sorting
0100                 if (cnt == b + 1)
0101                 {
0102                     std::sort(this->heights.begin(), this->heights.end());
0103                 }
0104             }
0105             else
0106             {
0107                 // find cell k such that heights[k-1] <= args[sample] < heights[k] and adjust extreme values
0108                 if (args[sample] < this->heights[0])
0109                 {
0110                     this->heights[0] = args[sample];
0111                     sample_cell = 1;
0112                 }
0113                 else if (this->heights[b] <= args[sample])
0114                 {
0115                     this->heights[b] = args[sample];
0116                     sample_cell = b;
0117                 }
0118                 else
0119                 {
0120                     typename array_type::iterator it;
0121                     it = std::upper_bound(
0122                         this->heights.begin()
0123                       , this->heights.end()
0124                       , args[sample]
0125                     );
0126
0127                     sample_cell = std::distance(this->heights.begin(), it);
0128                 }
0129
0130                 // increment positions of markers above sample_cell
0131                 for (std::size_t i = sample_cell; i < b + 1; ++i)
0132                 {
0133                     ++this->actual_positions[i];
0134                 }
0135
0136                 // update desired position of markers 2 to num_cells + 1
0137                 // (desired position of first marker is always 1)
0138                 for (std::size_t i = 1; i < b + 1; ++i)
0139                 {
0140                     this->desired_positions[i] += this->positions_increments[i];
0141                 }
0142
0143                 // adjust heights of markers 2 to num_cells if necessary
0144                 for (std::size_t i = 1; i < b; ++i)
0145                 {
0146                     // offset to desire position
0147                     float_type d = this->desired_positions[i] - this->actual_positions[i];
0148
0149                     // offset to next position
0150                     float_type dp = this->actual_positions[i + 1] - this->actual_positions[i];
0151
0152                     // offset to previous position
0153                     float_type dm = this->actual_positions[i - 1] - this->actual_positions[i];
0154
0155                     // height ds
0156                     float_type hp = (this->heights[i + 1] - this->heights[i]) / dp;
0157                     float_type hm = (this->heights[i - 1] - this->heights[i]) / dm;
0158
0159                     if ( ( d >= 1. && dp > 1. ) || ( d <= -1. && dm < -1. ) )
0160                     {
0161                         short sign_d = static_cast<short>(d / std::abs(d));
0162
0163                         // try adjusting heights[i] using p-squared formula
0164                         float_type h = this->heights[i] + sign_d / (dp - dm) * ( (sign_d - dm) * hp + (dp - sign_d) * hm );
0165
0166                         if ( this->heights[i - 1] < h && h < this->heights[i + 1] )
0167                         {
0168                             this->heights[i] = h;
0169                         }
0170                         else
0171                         {
0172                             // use linear formula
0173                             if (d>0)
0174                             {
0175                                 this->heights[i] += hp;
0176                             }
0177                             if (d<0)
0178                             {
0179                                 this->heights[i] -= hm;
0180                             }
0181                         }
0182                         this->actual_positions[i] += sign_d;
0183                     }
0184                 }
0185             }
0186         }
0187
0188         template<typename Args>
0189         result_type result(Args const &args) const
0190         {
0191             if (this->is_dirty)
0192             {
0193                 this->is_dirty = false;
0194
0195                 // creates a vector of std::pair where each pair i holds
0196                 // the values heights[i] (x-axis of histogram) and
0197                 // actual_positions[i] / cnt (y-axis of histogram)
0198
0199                 std::size_t cnt = count(args);
0200
0201                 for (std::size_t i = 0; i < this->histogram.size(); ++i)
0202                 {
0203                     this->histogram[i] = std::make_pair(this->heights[i], numeric::fdiv(this->actual_positions[i], cnt));
0204                 }
0205             }
0206             //return histogram;
0207             return make_iterator_range(this->histogram);
0208         }
0209
0210         // make this accumulator serializeable
0211         // TODO split to save/load and check on parameters provided in ctor
0212         template<class Archive>
0213         void serialize(Archive & ar, const unsigned int file_version)
0214         {
0215             ar & num_cells;
0216             ar & heights;
0217             ar & actual_positions;
0218             ar & desired_positions;
0219             ar & positions_increments;
0220             ar & histogram;
0221             ar & is_dirty;
0222         }
0223
0224     private:
0225         std::size_t num_cells;            // number of cells b
0226         array_type  heights;              // q_i
0227         array_type  actual_positions;     // n_i
0228         array_type  desired_positions;    // n'_i
0229         array_type  positions_increments; // dn'_i
0230         mutable histogram_type histogram; // histogram
0231         mutable bool is_dirty;
0232     };
0233
0234 } // namespace detail
0235
0236 ///////////////////////////////////////////////////////////////////////////////
0237 // tag::p_square_cumulative_distribution
0238 //
0239 namespace tag
0240 {
0241     struct p_square_cumulative_distribution
0242       : depends_on<count>
0243       , p_square_cumulative_distribution_num_cells
0244     {
0245         /// INTERNAL ONLY
0246         ///
0247         typedef accumulators::impl::p_square_cumulative_distribution_impl<mpl::_1> impl;
0248     };
0249 }
0250
0251 ///////////////////////////////////////////////////////////////////////////////
0252 // extract::p_square_cumulative_distribution
0253 //
0254 namespace extract
0255 {
0256     extractor<tag::p_square_cumulative_distribution> const p_square_cumulative_distribution = {};
0257
0258     BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_cumulative_distribution)
0259 }
0260
0261 using extract::p_square_cumulative_distribution;
0262
0263 // So that p_square_cumulative_distribution can be automatically substituted with
0264 // weighted_p_square_cumulative_distribution when the weight parameter is non-void
0265 template<>
0266 struct as_weighted_feature<tag::p_square_cumulative_distribution>
0267 {
0268     typedef tag::weighted_p_square_cumulative_distribution type;
0269 };
0270
0271 template<>
0272 struct feature_of<tag::weighted_p_square_cumulative_distribution>
0273   : feature_of<tag::p_square_cumulative_distribution>
0274 {
0275 };
0276
0277 }} // namespace boost::accumulators
0278
0279 #endif