Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
ComputationImpl.icpp
Go to the documentation of this file.
1namespace Euclid {
2namespace Histogram {
3
4template <typename VarType, typename WeightType>
5template <typename BinType>
6void 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
20template <typename VarType, typename WeightType>
21template <typename BinType>
22std::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 */
73template <typename BinType, typename IterType>
74struct 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 */
87template <typename BinType, typename IterType>
88inline 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 */
95template <typename BinType, typename IterType>
96inline void computeBinsForwarder(BinType&, IterType, IterType, std::false_type) {}
97
98template <typename VarType, typename WeightType>
99template <typename BinType>
100template <typename IterType, typename WeightIterType>
101void 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