Alexandria  2.18
Please provide a description of the project.
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>
6 void Histogram<VarType, WeightType>::ComputationImpl<BinType>::clip(VarType min, VarType max) {
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 template <typename VarType, typename WeightType>
21 template <typename BinType>
22 std::tuple<VarType, VarType, VarType> Histogram<VarType, WeightType>::ComputationImpl<BinType>::getStats() const {
23  VarType total = 0, total_count = 0;
24  VarType sigma = 0;
25 
26  // Find the mean and standard deviation in one go
27  for (auto i = m_clip_left; i <= m_clip_right; ++i) {
28  auto center = m_binning.getBin(i);
29  total += (*m_counts)[i] * center;
30  total_count += (*m_counts)[i];
31  sigma += (*m_counts)[i] * center * center;
32  }
33 
34  VarType mean = total / total_count;
35  sigma = sigma / total_count - mean * mean;
36  if (sigma > 0)
37  sigma = std::sqrt(sigma);
38  else
39  sigma = 0;
40 
41  // Find the median
42  WeightType low_sum = 0., high_sum = 0.;
43  auto low_i = m_clip_left, high_i = m_clip_right;
44  while (low_i <= high_i) {
45  if (low_sum < high_sum) {
46  low_sum += (*m_counts)[low_i++];
47  } else {
48  high_sum += (*m_counts)[high_i--];
49  }
50  }
51 
52  assert(low_sum + high_sum == total_count);
53 
54  VarType median;
55  if (high_i >= 0) {
56  auto edges = m_binning.getBinEdges(high_i + 1);
57  auto bin_width = (edges.second - edges.first);
58  auto max_counts = std::max((*m_counts)[low_i], (*m_counts)[high_i]);
59  median = edges.first + bin_width * (high_sum - low_sum) / (2.0 * max_counts);
60  } else {
61  median = m_binning.getBin(0);
62  }
63 
64  return std::make_tuple(mean, median, sigma);
65 }
66 
67 /**
68  * This class as a constexpr static member "value" which will be 'true' iff BinType has a method
69  * computeBins that can receive two instances of IterType as parameters
70  * @tparam BinType
71  * @tparam IterType
72  */
73 template <typename BinType, typename IterType>
74 struct HasComputeBins {
75  template <typename U, typename = decltype(std::declval<U>().computeBins(std::declval<IterType>(), std::declval<IterType>()))>
76  struct SFINAE;
77  template <typename U>
78  static char Test(SFINAE<U>*);
79  template <typename U>
80  static int Test(...);
81  static constexpr bool value = sizeof(Test<BinType>(0)) == sizeof(char);
82 };
83 
84 /**
85  * This method is called if BinType has computeBins
86  */
87 template <typename BinType, typename IterType>
88 inline void computeBinsForwarder(BinType& binning, IterType begin, IterType end, std::true_type) {
89  binning.computeBins(begin, end);
90 }
91 
92 /**
93  * This method is called if BinType does not have computeBins
94  */
95 template <typename BinType, typename IterType>
96 inline void computeBinsForwarder(BinType&, IterType, IterType, std::false_type) {}
97 
98 template <typename VarType, typename WeightType>
99 template <typename BinType>
100 template <typename IterType, typename WeightIterType>
101 void Histogram<VarType, WeightType>::ComputationImpl<BinType>::computeBins(IterType begin, IterType end, WeightIterType wbegin) {
102  // This trick should allow the compiler to know the actual binning type, so if they
103  // override methods with *final*, we can skip indirections via vtable
104  computeBinsForwarder(m_binning, begin, end, std::integral_constant<bool, HasComputeBins<BinType, IterType>::value>());
105 
106  m_clip_left = 0;
107  m_clip_right = m_binning.getBinCount() - 1;
108  m_counts->resize(m_binning.getBinCount());
109 
110  ssize_t nbins = m_counts->size();
111  auto i = begin;
112  auto wi = wbegin;
113  for (; i != end; ++i, ++wi) {
114  auto bin = m_binning.getBinIndex(*i);
115  if (bin >= 0 && bin < nbins) {
116  (*m_counts)[bin] += *wi;
117  }
118  }
119 }
120 
121 } // end of namespace Histogram
122 } // end of namespace Euclid
Euclid
Definition: InstOrRefHolder.h:29
Euclid::Histogram::Histogram
Definition: Histogram.h:137