1 #pragma once
2 
3 /// \file Means.h
4 /// \brief Computing minimum, maximum and mean value of floats.
5 /// \author Pavel Sevecek (sevecek at sirrah.troja.mff.cuni.cz)
6 /// \date 2016-2021
7 
8 #include "objects/geometry/Generic.h"
9 #include "objects/wrappers/Interval.h"
10 
11 NAMESPACE_SPH_BEGIN
12 
13 /// \brief Generalized mean with fixed (compile-time) power.
14 template <int Power>
15 class GeneralizedMean {
16 private:
17     double sum = 0.; // using double to limit round-off errors in summing
18     Size weight = 0;
19 
20 public:
21     GeneralizedMean() = default;
22 
accumulate(const Float value)23     INLINE void accumulate(const Float value) {
24         sum += pow<Power>(value);
25         weight++;
26     }
27 
28     /// Removes all values from the set.
reset()29     INLINE void reset() {
30         sum = 0.;
31         weight = 0;
32     }
33 
compute()34     INLINE Float compute() const {
35         return Float(pow(sum / weight, 1. / Power));
36     }
37 
count()38     INLINE Size count() const {
39         return weight;
40     }
41 
42     friend std::ostream& operator<<(std::ostream& stream, const GeneralizedMean& stats) {
43         stream << stats.compute();
44         return stream;
45     }
46 };
47 
48 /// Geometric mean has to be specialized.
49 template <>
50 class GeneralizedMean<0> {
51 private:
52     double sum = 1.;
53     Size weight = 0;
54 
55 public:
56     GeneralizedMean() = default;
57 
accumulate(const Float value)58     INLINE void accumulate(const Float value) {
59         sum *= value;
60         weight++;
61     }
62 
reset()63     INLINE void reset() {
64         sum = 1.;
65         weight = 0;
66     }
67 
compute()68     INLINE Float compute() const {
69         return Float(pow(sum, 1. / weight));
70     }
71 
count()72     INLINE Size count() const {
73         return weight;
74     }
75 
76     friend std::ostream& operator<<(std::ostream& stream, const GeneralizedMean& stats) {
77         stream << stats.compute();
78         return stream;
79     }
80 };
81 
82 
83 /// Aliases
84 using ArithmeticMean = GeneralizedMean<1>;
85 using HarmonicMean = GeneralizedMean<-1>;
86 using GeometricMean = GeneralizedMean<0>;
87 
88 
89 /// \brief Generalized mean with positive (runtime) power.
90 ///
91 /// Cannot be used to compute geometric mean. If constructed with non-positive power, it issues an assert.
92 class PositiveMean {
93 protected:
94     double sum = 0.;
95     Size weight = 0;
96     Float power;
97 
98 public:
PositiveMean(const Float power)99     PositiveMean(const Float power)
100         : power(power) {
101         SPH_ASSERT(power > 0._f);
102     }
103 
accumulate(const Float value)104     INLINE void accumulate(const Float value) {
105         sum += powFastest(value, power);
106         weight++;
107     }
108 
accumulate(const PositiveMean & other)109     INLINE void accumulate(const PositiveMean& other) {
110         SPH_ASSERT(power == other.power); // it only makes sense to sum up same means
111         sum += other.sum;
112         weight += other.weight;
113     }
114 
reset()115     INLINE void reset() {
116         sum = 0.;
117         weight = 0;
118     }
119 
compute()120     INLINE Float compute() const {
121         return Float(powFastest(Float(sum / weight), 1._f / power));
122     }
123 
count()124     INLINE Size count() const {
125         return weight;
126     }
127 
128     friend std::ostream& operator<<(std::ostream& stream, const PositiveMean& stats) {
129         stream << stats.compute();
130         return stream;
131     }
132 };
133 
134 /// \brief Generalized mean with negative (runtime) power.
135 ///
136 /// Cannot be used to compute geometric mean. If constructed with non-negative power, it issues an assert.
137 class NegativeMean : public PositiveMean {
138 public:
NegativeMean(const Float power)139     NegativeMean(const Float power)
140         : PositiveMean(-power) {}
141 
accumulate(const Float value)142     INLINE void accumulate(const Float value) {
143         SPH_ASSERT(value > 0._f, value);
144         const Float p = pow(value, power);
145         if (p == INFTY) {
146             weight++; // just increase weight
147         } else if (p > 0._f) {
148             sum += 1._f / p;
149             weight++;
150         }
151     }
152 
accumulate(const NegativeMean & other)153     INLINE void accumulate(const NegativeMean& other) {
154         SPH_ASSERT(power == other.power);
155         sum += other.sum;
156         weight += other.weight;
157     }
158 
compute()159     INLINE Float compute() const {
160         Float avg = Float(sum / weight);
161         SPH_ASSERT(isReal(avg), avg, sum, weight);
162         Float avgPow = pow(avg, 1._f / power);
163         if (avgPow == 0._f) {
164             return INFINITY;
165         } else {
166             return 1._f / avgPow;
167         }
168     }
169 };
170 
171 /// Helper class for statistics, accumulating minimal, maximal and mean value of a set of numbers.
172 class MinMaxMean {
173 private:
174     Interval minMax;
175     ArithmeticMean avg;
176 
177 public:
178     MinMaxMean() = default;
179 
accumulate(const Float value)180     INLINE void accumulate(const Float value) {
181         avg.accumulate(value);
182         minMax.extend(value);
183     }
184 
185     /// Removes all values from the set.
reset()186     INLINE void reset() {
187         avg.reset();
188         minMax = Interval();
189     }
190 
mean()191     INLINE Float mean() const {
192         return avg.compute();
193     }
194 
min()195     INLINE Float min() const {
196         return minMax.lower();
197     }
198 
max()199     INLINE Float max() const {
200         return minMax.upper();
201     }
202 
range()203     INLINE Interval range() const {
204         return minMax;
205     }
206 
count()207     INLINE Size count() const {
208         return avg.count();
209     }
210 
211     template <typename TStream>
212     friend TStream& operator<<(TStream& stream, const MinMaxMean& stats) {
213         stream << "average = " << stats.mean() << "  (min = " << stats.min() << ", max = " << stats.max()
214                << ")";
215         return stream;
216     }
217 };
218 
219 NAMESPACE_SPH_END
220