/Src/Dependencies/Boost/boost/accumulators/statistics/p_square_cumulative_distribution.hpp

  1///////////////////////////////////////////////////////////////////////////////
  2// p_square_cumulative_distribution.hpp
  3//
  4//  Copyright 2005 Daniel Egloff, Olivier Gygi. Distributed under the Boost
  5//  Software License, Version 1.0. (See accompanying file
  6//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  7
  8#ifndef BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_CUMULATIVE_DISTRIBUTION_HPP_DE_01_01_2006
  9#define BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_CUMULATIVE_DISTRIBUTION_HPP_DE_01_01_2006
 10
 11#include <vector>
 12#include <functional>
 13#include <boost/parameter/keyword.hpp>
 14#include <boost/range.hpp>
 15#include <boost/mpl/placeholders.hpp>
 16#include <boost/accumulators/framework/accumulator_base.hpp>
 17#include <boost/accumulators/framework/extractor.hpp>
 18#include <boost/accumulators/numeric/functional.hpp>
 19#include <boost/accumulators/framework/parameters/sample.hpp>
 20#include <boost/accumulators/statistics_fwd.hpp>
 21#include <boost/accumulators/statistics/count.hpp>
 22
 23namespace boost { namespace accumulators
 24{
 25///////////////////////////////////////////////////////////////////////////////
 26// num_cells named parameter
 27//
 28BOOST_PARAMETER_NESTED_KEYWORD(tag, p_square_cumulative_distribution_num_cells, num_cells)
 29
 30namespace impl
 31{
 32    ///////////////////////////////////////////////////////////////////////////////
 33    // p_square_cumulative_distribution_impl
 34    //  cumulative_distribution calculation (as histogram)
 35    /**
 36        @brief Histogram calculation of the cumulative distribution with the \f$P^2\f$ algorithm
 37
 38        A histogram of the sample cumulative distribution is computed dynamically without storing samples
 39        based on the \f$ P^2 \f$ algorithm. The returned histogram has a specifiable amount (num_cells)
 40        equiprobable (and not equal-sized) cells.
 41
 42        For further details, see
 43
 44        R. Jain and I. Chlamtac, The P^2 algorithmus for dynamic calculation of quantiles and
 45        histograms without storing observations, Communications of the ACM,
 46        Volume 28 (October), Number 10, 1985, p. 1076-1085.
 47
 48        @param p_square_cumulative_distribution_num_cells.
 49    */
 50    template<typename Sample>
 51    struct p_square_cumulative_distribution_impl
 52      : accumulator_base
 53    {
 54        typedef typename numeric::functional::average<Sample, std::size_t>::result_type float_type;
 55        typedef std::vector<float_type> array_type;
 56        typedef std::vector<std::pair<float_type, float_type> > histogram_type;
 57        // for boost::result_of
 58        typedef iterator_range<typename histogram_type::iterator> result_type;
 59
 60        template<typename Args>
 61        p_square_cumulative_distribution_impl(Args const &args)
 62          : num_cells(args[p_square_cumulative_distribution_num_cells])
 63          , heights(num_cells + 1)
 64          , actual_positions(num_cells + 1)
 65          , desired_positions(num_cells + 1)
 66          , positions_increments(num_cells + 1)
 67          , histogram(num_cells + 1)
 68          , is_dirty(true)
 69        {
 70            std::size_t b = this->num_cells;
 71
 72            for (std::size_t i = 0; i < b + 1; ++i)
 73            {
 74                this->actual_positions[i] = i + 1.;
 75                this->desired_positions[i] = i + 1.;
 76                this->positions_increments[i] = numeric::average(i, b);
 77            }
 78        }
 79
 80        template<typename Args>
 81        void operator ()(Args const &args)
 82        {
 83            this->is_dirty = true;
 84
 85            std::size_t cnt = count(args);
 86            std::size_t sample_cell = 1; // k
 87            std::size_t b = this->num_cells;
 88
 89            // accumulate num_cells + 1 first samples
 90            if (cnt <= b + 1)
 91            {
 92                this->heights[cnt - 1] = args[sample];
 93
 94                // complete the initialization of heights by sorting
 95                if (cnt == b + 1)
 96                {
 97                    std::sort(this->heights.begin(), this->heights.end());
 98                }
 99            }
100            else
101            {
102                // find cell k such that heights[k-1] <= args[sample] < heights[k] and adjust extreme values
103                if (args[sample] < this->heights[0])
104                {
105                    this->heights[0] = args[sample];
106                    sample_cell = 1;
107                }
108                else if (this->heights[b] <= args[sample])
109                {
110                    this->heights[b] = args[sample];
111                    sample_cell = b;
112                }
113                else
114                {
115                    typename array_type::iterator it;
116                    it = std::upper_bound(
117                        this->heights.begin()
118                      , this->heights.end()
119                      , args[sample]
120                    );
121
122                    sample_cell = std::distance(this->heights.begin(), it);
123                }
124
125                // increment positions of markers above sample_cell
126                for (std::size_t i = sample_cell; i < b + 1; ++i)
127                {
128                    ++this->actual_positions[i];
129                }
130
131                // update desired position of markers 2 to num_cells + 1
132                // (desired position of first marker is always 1)
133                for (std::size_t i = 1; i < b + 1; ++i)
134                {
135                    this->desired_positions[i] += this->positions_increments[i];
136                }
137
138                // adjust heights of markers 2 to num_cells if necessary
139                for (std::size_t i = 1; i < b; ++i)
140                {
141                    // offset to desire position
142                    float_type d = this->desired_positions[i] - this->actual_positions[i];
143
144                    // offset to next position
145                    float_type dp = this->actual_positions[i + 1] - this->actual_positions[i];
146
147                    // offset to previous position
148                    float_type dm = this->actual_positions[i - 1] - this->actual_positions[i];
149
150                    // height ds
151                    float_type hp = (this->heights[i + 1] - this->heights[i]) / dp;
152                    float_type hm = (this->heights[i - 1] - this->heights[i]) / dm;
153
154                    if ( ( d >= 1. && dp > 1. ) || ( d <= -1. && dm < -1. ) )
155                    {
156                        short sign_d = static_cast<short>(d / std::abs(d));
157
158                        // try adjusting heights[i] using p-squared formula
159                        float_type h = this->heights[i] + sign_d / (dp - dm) * ( (sign_d - dm) * hp + (dp - sign_d) * hm );
160
161                        if ( this->heights[i - 1] < h && h < this->heights[i + 1] )
162                        {
163                            this->heights[i] = h;
164                        }
165                        else
166                        {
167                            // use linear formula
168                            if (d>0)
169                            {
170                                this->heights[i] += hp;
171                            }
172                            if (d<0)
173                            {
174                                this->heights[i] -= hm;
175                            }
176                        }
177                        this->actual_positions[i] += sign_d;
178                    }
179                }
180            }
181        }
182
183        template<typename Args>
184        result_type result(Args const &args) const
185        {
186            if (this->is_dirty)
187            {
188                this->is_dirty = false;
189
190                // creates a vector of std::pair where each pair i holds
191                // the values heights[i] (x-axis of histogram) and
192                // actual_positions[i] / cnt (y-axis of histogram)
193
194                std::size_t cnt = count(args);
195
196                for (std::size_t i = 0; i < this->histogram.size(); ++i)
197                {
198                    this->histogram[i] = std::make_pair(this->heights[i], numeric::average(this->actual_positions[i], cnt));
199                }
200            }
201            //return histogram;
202            return make_iterator_range(this->histogram);
203        }
204
205    private:
206        std::size_t num_cells;            // number of cells b
207        array_type  heights;              // q_i
208        array_type  actual_positions;     // n_i
209        array_type  desired_positions;    // n'_i
210        array_type  positions_increments; // dn'_i
211        mutable histogram_type histogram; // histogram
212        mutable bool is_dirty;
213    };
214
215} // namespace detail
216
217///////////////////////////////////////////////////////////////////////////////
218// tag::p_square_cumulative_distribution
219//
220namespace tag
221{
222    struct p_square_cumulative_distribution
223      : depends_on<count>
224      , p_square_cumulative_distribution_num_cells
225    {
226        /// INTERNAL ONLY
227        ///
228        typedef accumulators::impl::p_square_cumulative_distribution_impl<mpl::_1> impl;
229    };
230}
231
232///////////////////////////////////////////////////////////////////////////////
233// extract::p_square_cumulative_distribution
234//
235namespace extract
236{
237    extractor<tag::p_square_cumulative_distribution> const p_square_cumulative_distribution = {};
238
239    BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_cumulative_distribution)
240}
241
242using extract::p_square_cumulative_distribution;
243
244// So that p_square_cumulative_distribution can be automatically substituted with
245// weighted_p_square_cumulative_distribution when the weight parameter is non-void
246template<>
247struct as_weighted_feature<tag::p_square_cumulative_distribution>
248{
249    typedef tag::weighted_p_square_cumulative_distribution type;
250};
251
252template<>
253struct feature_of<tag::weighted_p_square_cumulative_distribution>
254  : feature_of<tag::p_square_cumulative_distribution>
255{
256};
257
258}} // namespace boost::accumulators
259
260#endif