yat
0.8.3pre
|
00001 #ifndef _theplu_yat_statistics_percentiler_ 00002 #define _theplu_yat_statistics_percentiler_ 00003 00004 // $Id: Percentiler.h 2470 2011-04-12 03:21:45Z peter $ 00005 00006 /* 00007 Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson 00008 Copyright (C) 2010, 2011 Peter Johansson 00009 00010 This file is part of the yat library, http://dev.thep.lu.se/yat 00011 00012 The yat library is free software; you can redistribute it and/or 00013 modify it under the terms of the GNU General Public License as 00014 published by the Free Software Foundation; either version 3 of the 00015 License, or (at your option) any later version. 00016 00017 The yat library is distributed in the hope that it will be useful, 00018 but WITHOUT ANY WARRANTY; without even the implied warranty of 00019 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00020 General Public License for more details. 00021 00022 You should have received a copy of the GNU General Public License 00023 along with yat. If not, see <http://www.gnu.org/licenses/>. 00024 */ 00025 00026 #include "yat/utility/concept_check.h" 00027 #include "yat/utility/DataWeight.h" 00028 #include "yat/utility/iterator_traits.h" 00029 #include "yat/utility/yat_assert.h" 00030 #include "yat/utility/WeightIterator.h" 00031 00032 #include <boost/concept_check.hpp> 00033 00034 #include <algorithm> 00035 #include <cmath> 00036 #include <limits> 00037 #include <numeric> 00038 #include <stdexcept> 00039 #include <vector> 00040 00041 namespace theplu { 00042 namespace yat { 00043 namespace statistics { 00044 00050 class Percentiler 00051 { 00052 public: 00060 Percentiler(double perc=50, bool sorted=false); 00061 00096 template<typename RandomAccessIterator> 00097 double operator()(RandomAccessIterator first, 00098 RandomAccessIterator last) const 00099 { 00100 BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>)); 00101 BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>)); 00102 if (first==last) 00103 return std::numeric_limits<double>::quiet_NaN(); 00104 return calculate(first, last, sorted_, 00105 typename utility::weighted_iterator_traits<RandomAccessIterator>::type()); 00106 } 00107 private: 00108 double perc_; 00109 bool sorted_; 00110 00111 // unweighted version 00112 template<typename RandomAccessIterator> 00113 double calculate(RandomAccessIterator first, RandomAccessIterator last, 00114 bool sorted, utility::unweighted_iterator_tag tag) const; 00115 00116 // weighted version 00117 template<typename RandomAccessIterator> 00118 double calculate(RandomAccessIterator first, RandomAccessIterator last, 00119 bool sorted, utility::weighted_iterator_tag tag) const; 00120 00121 // using compiler generated copy 00122 //Percentiler(const Percentiler&); 00123 //Percentiler& operator=(const Percentiler&); 00124 00125 }; 00126 00127 // template implementations 00128 // ///////////////////////// 00129 00130 // unweighted version 00131 template<typename RandomAccessIterator> 00132 double 00133 Percentiler::calculate(RandomAccessIterator first, 00134 RandomAccessIterator last, 00135 bool sorted, 00136 utility::unweighted_iterator_tag tag) const 00137 { 00138 // range is one value only is a special case 00139 if (first+1 == last) 00140 *first; 00141 if (sorted) { 00142 size_t n = last - first; 00143 // have to take care of this special case 00144 if (perc_>= 100.0) 00145 return *(--last); 00146 if (perc_<= 0.0) 00147 return *first; 00148 double j = n * perc_ / 100.0; 00149 size_t i = static_cast<size_t>(j); 00150 if (i==j) 00151 return (first[i]+first[i-1])/2; 00152 return first[i]; 00153 } 00154 00155 std::vector<double> v_copy; 00156 v_copy.reserve(std::distance(first,last)); 00157 std::copy(first, last, std::back_inserter(v_copy)); 00158 std::sort(v_copy.begin(), v_copy.end()); 00159 return calculate(v_copy.begin(), v_copy.end(), true, tag); 00160 } 00161 00162 00163 // weighted version 00164 template<typename RandomAccessIterator> 00165 double Percentiler::calculate(RandomAccessIterator first, 00166 RandomAccessIterator last, 00167 bool sorted, 00168 utility::weighted_iterator_tag tag) const 00169 { 00170 if (sorted) { 00171 utility::iterator_traits<RandomAccessIterator> trait; 00172 std::vector<double> accum_w; 00173 accum_w.reserve(last-first); 00174 std::partial_sum(weight_iterator(first), 00175 weight_iterator(last), 00176 std::back_inserter(accum_w)); 00177 00178 double w_bound=perc_/100.0*accum_w.back(); 00179 std::vector<double>::const_iterator upper(accum_w.begin()); 00180 double margin=1e-10; 00181 while (upper!=accum_w.end() && *upper <= w_bound+margin) 00182 ++upper; 00183 while (upper!=accum_w.begin() && 00184 (upper==accum_w.end() || 00185 trait.weight(first+(upper-accum_w.begin()))==0.0)) 00186 --upper; 00187 std::vector<double>::const_iterator lower(upper); 00188 while ( *(lower-1)>=w_bound-margin && lower>accum_w.begin()) 00189 --lower; 00190 00191 return (trait.data(first+(upper-accum_w.begin()))+ 00192 trait.data(first+(lower-accum_w.begin())))/2; 00193 } 00194 00195 std::vector<utility::DataWeight> v_copy; 00196 v_copy.reserve(last-first); 00197 std::copy(first, last, std::back_inserter(v_copy)); 00198 std::sort(v_copy.begin(), v_copy.end()); 00199 return calculate(v_copy.begin(), v_copy.end(), true, tag); 00200 } 00201 00202 00203 }}} // of namespace statistics, yat, and theplu 00204 00205 #endif