00001 #ifndef _theplu_yat_statistics_percentiler_
00002 #define _theplu_yat_statistics_percentiler_
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "yat/utility/DataWeight.h"
00027 #include "yat/utility/iterator_traits.h"
00028 #include "yat/utility/yat_assert.h"
00029 #include "yat/utility/WeightIterator.h"
00030
00031 #include <algorithm>
00032 #include <cassert>
00033 #include <cmath>
00034 #include <numeric>
00035 #include <stdexcept>
00036 #include <vector>
00037
00038 #include <iostream>
00039
00040 namespace theplu {
00041 namespace yat {
00042 namespace statistics {
00043
00049 class Percentiler
00050 {
00051 public:
00059 Percentiler(double perc=50, bool sorted=false);
00060
00092 template<typename RandomAccessIterator>
00093 double operator()(RandomAccessIterator first,
00094 RandomAccessIterator last) const
00095 {
00096 return calculate(first, last, sorted_,
00097 typename utility::weighted_iterator_traits<RandomAccessIterator>::type());
00098 }
00099 private:
00100 double perc_;
00101 bool sorted_;
00102
00103
00104 template<typename RandomAccessIterator>
00105 double calculate(RandomAccessIterator first, RandomAccessIterator last,
00106 bool sorted, utility::unweighted_iterator_tag tag) const;
00107
00108
00109 template<typename RandomAccessIterator>
00110 double calculate(RandomAccessIterator first, RandomAccessIterator last,
00111 bool sorted, utility::weighted_iterator_tag tag) const;
00112
00113
00114
00115
00116
00117 };
00118
00119
00120
00121
00122
00123 template<typename RandomAccessIterator>
00124 double
00125 Percentiler::calculate(RandomAccessIterator first,
00126 RandomAccessIterator last,
00127 bool sorted,
00128 utility::unweighted_iterator_tag tag) const
00129 {
00130
00131 if (first+1 == last)
00132 *first;
00133 if (sorted) {
00134 size_t n = last - first;
00135
00136 if (perc_>= 100.0)
00137 return *(--last);
00138 if (perc_<= 0.0)
00139 return *first;
00140 double j = n * perc_ / 100.0;
00141 size_t i = static_cast<size_t>(j);
00142 if (i==j)
00143 return (first[i]+first[i-1])/2;
00144 return first[i];
00145 }
00146
00147 std::vector<double> v_copy;
00148 v_copy.reserve(std::distance(first,last));
00149 std::copy(first, last, std::back_inserter(v_copy));
00150 std::sort(v_copy.begin(), v_copy.end());
00151 return calculate(v_copy.begin(), v_copy.end(), true, tag);
00152 }
00153
00154
00155
00156 template<typename RandomAccessIterator>
00157 double Percentiler::calculate(RandomAccessIterator first,
00158 RandomAccessIterator last,
00159 bool sorted,
00160 utility::weighted_iterator_tag tag) const
00161 {
00162 if (sorted) {
00163 utility::iterator_traits<RandomAccessIterator> trait;
00164 std::vector<double> accum_w;
00165 accum_w.reserve(last-first);
00166 std::partial_sum(weight_iterator(first),
00167 weight_iterator(last),
00168 std::back_inserter(accum_w));
00169
00170 double w_bound=perc_/100.0*accum_w.back();
00171 std::vector<double>::const_iterator upper(accum_w.begin());
00172 double margin=1e-10;
00173 while (*upper <= w_bound+margin && upper!=accum_w.end())
00174 ++upper;
00175 if (upper==accum_w.end())
00176 --upper;
00177 std::vector<double>::const_iterator lower(upper);
00178 while ( *(lower-1)>=w_bound-margin && lower>accum_w.begin())
00179 --lower;
00180
00181 return (trait.data(first+(upper-accum_w.begin()))+
00182 trait.data(first+(lower-accum_w.begin())))/2;
00183 }
00184
00185 std::vector<utility::DataWeight> v_copy;
00186 v_copy.reserve(last-first);
00187 utility::iterator_traits<RandomAccessIterator> traits;
00188 for ( ; first!=last; ++first)
00189 if (traits.weight(first))
00190 v_copy.push_back(*first);
00191 std::sort(v_copy.begin(), v_copy.end());
00192 return calculate(v_copy.begin(), v_copy.end(), true, tag);
00193 }
00194
00195
00196 }}}
00197
00198 #endif