Statistics/interface/hsm_1d.h

0001 #ifndef CommonTools_Statistics_hsm_1d_h
0002 #define CommonTools_Statistics_hsm_1d_h
0003
0004 #include "CommonTools/Statistics/interface/StatisticsException.h"
0005
0006 #include <vector>
0007 #include <cmath>
0008 #include <algorithm>
0009
0010 template <class T>
0011 std::pair<typename std::vector<T>::iterator, typename std::vector<T>::iterator> hsm_1d_returning_iterators(
0012     std::vector<T>& values) {
0013   sort(values.begin(), values.end());
0014
0015   typename std::vector<T>::iterator begin = values.begin();
0016   typename std::vector<T>::iterator end = values.end() - 1;
0017
0018   while (true) {
0019     const int size = (int)(end - begin);
0020     // const int stepsize = size / 2;
0021     int stepsize = static_cast<int>(floor(.4999999999 + size * .5));
0022     if (stepsize == 0)
0023       stepsize = 1;
0024
0025     end = begin + stepsize;
0026     T div = (T)fabs((T)(*end) - (T)(*begin));
0027     typename std::vector<T>::iterator new_begin = begin;
0028     typename std::vector<T>::iterator new_end = end;
0029
0030     for (typename std::vector<T>::iterator i = (begin + 1); i != (begin + size - stepsize + 1); ++i) {
0031       if (fabs((*i) - (*(i + stepsize))) < div) {
0032         new_begin = i;
0033         new_end = i + stepsize;
0034         div = (T)((T)(*new_end) - (T)(*new_begin));
0035       };
0036     };
0037     begin = new_begin;
0038     end = new_end;
0039     if (size < 4)
0040       break;
0041   };
0042
0043   std::pair<typename std::vector<T>::iterator, typename std::vector<T>::iterator> ret(begin, end);
0044   return ret;
0045 }
0046
0047 /**
0048  *  stepsize Sample Mode one dimension.
0049  *  T must be sortable, 'addable', and dividable by 2.
0050  */
0051
0052 template <class T>
0053 T hsm_1d(std::vector<T> values) {
0054   switch (values.size()) {
0055     case 0:
0056       throw StatisticsException("Hsm of empty vector undefined");
0057     case 1:
0058       return (T) * (values.begin());
0059     case 2:
0060       return (T)((*(values.begin()) + *(values.end() - 1)) / 2);
0061   };
0062   std::pair<typename std::vector<T>::iterator, typename std::vector<T>::iterator> itors =
0063       hsm_1d_returning_iterators(values);
0064   T ret = ((T)(*(itors.first)) + (T)(*(itors.second))) / (T)2.;
0065   return ret;
0066 }
0067
0068 #endif