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

http://hadesmem.googlecode.com/ · C++ Header · 260 lines · 171 code · 34 blank · 55 comment · 20 complexity · ad0a351999d1e8e0676d3fbbe6fb92bb MD5 · raw file

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