Alexandria  2.16
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ComputationImpl.icpp
Go to the documentation of this file.
1 namespace Euclid {
2 namespace Histogram {
3 
4 template<typename VarType, typename WeightType>
5 template<typename BinType>
7 
8  if (min > max)
9  throw Elements::Exception("Clipping with min > max can not be done");
10 
11  auto min_bin = m_binning.getBinIndex(min);
12  if (min_bin > m_clip_left && min_bin <= m_clip_right)
13  m_clip_left = min_bin;
14 
15  auto max_bin = m_binning.getBinIndex(max);
16  if (max_bin >= m_clip_left && max_bin < m_clip_right)
17  m_clip_right = max_bin;
18 
19 }
20 
21 template<typename VarType, typename WeightType>
22 template<typename BinType>
24  VarType total = 0, total_count = 0;
25  VarType sigma = 0;
26 
27  // Find the mean and standard deviation in one go
28  for (auto i = m_clip_left; i <= m_clip_right; ++i) {
29  auto center = m_binning.getBin(i);
30  total += (*m_counts)[i] * center;
31  total_count += (*m_counts)[i];
32  sigma += (*m_counts)[i] * center * center;
33  }
34 
35  VarType mean = total / total_count;
36  sigma = sigma / total_count - mean * mean;
37  if (sigma > 0)
38  sigma = std::sqrt(sigma);
39  else
40  sigma = 0;
41 
42  // Find the median
43  WeightType low_sum = 0., high_sum = 0.;
44  auto low_i = m_clip_left, high_i = m_clip_right;
45  while (low_i <= high_i) {
46  if (low_sum < high_sum) {
47  low_sum += (*m_counts)[low_i++];
48  }
49  else {
50  high_sum += (*m_counts)[high_i--];
51  }
52  }
53 
54  assert(low_sum + high_sum == total_count);
55 
56  VarType median;
57  if (high_i >= 0) {
58  auto edges = m_binning.getBinEdges(high_i + 1);
59  auto bin_width = (edges.second - edges.first);
60  auto max_counts = std::max((*m_counts)[low_i], (*m_counts)[high_i]);
61  median = edges.first + bin_width * (high_sum - low_sum) / (2.0 * max_counts);
62  }
63  else {
64  median = m_binning.getBin(0);
65  }
66 
67  return std::make_tuple(mean, median, sigma);
68 }
69 
76 template<typename BinType, typename IterType>
78 {
79  template <typename U, typename = decltype(std::declval<U>().computeBins(std::declval<IterType>(), std::declval<IterType>()))> struct SFINAE;
80  template<typename U> static char Test(SFINAE<U>*);
81  template<typename U> static int Test(...);
82  static constexpr bool value = sizeof(Test<BinType>(0)) == sizeof(char);
83 };
84 
88 template<typename BinType, typename IterType>
89 inline void computeBinsForwarder(BinType& binning, IterType begin, IterType end, std::true_type) {
90  binning.computeBins(begin, end);
91 }
92 
96 template<typename BinType, typename IterType>
97 inline void computeBinsForwarder(BinType&, IterType, IterType, std::false_type) {
98 }
99 
100 template<typename VarType, typename WeightType>
101 template<typename BinType>
102 template<typename IterType, typename WeightIterType>
103 void Histogram<VarType, WeightType>::ComputationImpl<BinType>::computeBins(IterType begin, IterType end, WeightIterType wbegin) {
104  // This trick should allow the compiler to know the actual binning type, so if they
105  // override methods with *final*, we can skip indirections via vtable
107 
108  m_clip_left = 0;
109  m_clip_right = m_binning.getBinCount() - 1;
110  m_counts->resize(m_binning.getBinCount());
111 
112  ssize_t nbins = m_counts->size();
113  auto i = begin;
114  auto wi = wbegin;
115  for (; i != end; ++i, ++wi) {
116  auto bin = m_binning.getBinIndex(*i);
117  if (bin >= 0 && bin < nbins) {
118  (*m_counts)[bin] += *wi;
119  }
120  }
121 }
122 
123 } // end of namespace Histogram
124 } // end of namespace Euclid
void computeBinsForwarder(BinType &binning, IterType begin, IterType end, std::true_type)
T make_tuple(T...args)
static char Test(SFINAE< U > *)
T max(T...args)
T sqrt(T...args)