Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
Util.h
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Author: L. Moneta Tue Nov 14 15:44:38 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Utility functions for all ROOT Math classes
12
13#ifndef ROOT_Math_Util
14#define ROOT_Math_Util
15
16#include <algorithm>
17#include <chrono>
18#include <cmath>
19#include <functional>
20#include <limits>
21#include <numeric>
22#include <sstream>
23#include <string>
24
25
26// This can be protected against by defining ROOT_Math_VecTypes
27// This is only used for the R__HAS_VECCORE define
28// and a single VecCore function in EvalLog
29#ifndef ROOT_Math_VecTypes
30#include "Types.h"
31#endif
32
33
34// for defining unused variables in the interfaces
35// and have still them in the documentation
36#define MATH_UNUSED(var) (void)var
37
38
39namespace ROOT {
40
41namespace Math {
42
43/**
44 namespace defining Utility functions needed by mathcore
45*/
46namespace Util {
47
49
50public:
51 TimingScope(std::function<void(std::string const&)> printer, std::string const &message);
52
54
55private:
56 std::chrono::steady_clock::time_point fBegin;
57 std::function<void(std::string const&)> fPrinter;
58 const std::string fMessage;
59};
60
61 /**
62 Utility function for conversion to strings
63 */
64 template <class T>
65 std::string ToString(const T &val)
66 {
67 std::ostringstream buf;
68 buf << val;
69
70 std::string ret = buf.str();
71 return ret;
72 }
73
74 /// safe evaluation of log(x) with a protections against negative or zero argument to the log
75 /// smooth linear extrapolation below function values smaller than epsilon
76 /// (better than a simple cut-off)
77
78 template<class T>
79 inline T EvalLog(T x) {
80 static const T epsilon = T(2.0 * std::numeric_limits<double>::min());
81#ifdef R__HAS_VECCORE
82 T logval = vecCore::Blend<T>(x <= epsilon, x / epsilon + std::log(epsilon) - T(1.0), std::log(x));
83#else
84 T logval = x <= epsilon ? x / epsilon + std::log(epsilon) - T(1.0) : std::log(x);
85#endif
86 return logval;
87 }
88
89 } // end namespace Util
90
91 /// \class KahanSum
92 /// The Kahan summation is a compensated summation algorithm, which significantly reduces numerical errors
93 /// when adding a sequence of finite-precision floating point numbers.
94 /// This is done by keeping a separate running compensation (a variable to accumulate small errors).
95 ///
96 /// ### Auto-vectorisable accumulation
97 /// This class can internally use multiple accumulators (template parameter `N`).
98 /// When filled from a collection that supports index access from a *contiguous* block of memory,
99 /// compilers such as gcc, clang and icc can auto-vectorise the accumulation. This happens by cycling
100 /// through the internal accumulators based on the value of "`index % N`", so `N` accumulators can be filled from a block
101 /// of `N` numbers in a single instruction.
102 ///
103 /// The usage of multiple accumulators might slightly increase the precision in comparison to the single-accumulator version
104 /// with `N = 1`.
105 /// This depends on the order and magnitude of the numbers being accumulated. Therefore, in rare cases, the accumulation
106 /// result can change *in dependence of N*, even when the data are identical.
107 /// The magnitude of such differences is well below the precision of the floating point type, and will therefore mostly show
108 /// in the compensation sum(see Carry()). Increasing the number of accumulators therefore only makes sense to
109 /// speed up the accumulation, but not to increase precision.
110 ///
111 /// \param T The type of the values to be accumulated.
112 /// \param N Number of accumulators. Defaults to 1. Ideal values are the widths of a vector register on the relevant architecture.
113 /// Depending on the instruction set, good values are:
114 /// - AVX2-float: 8
115 /// - AVX2-double: 4
116 /// - AVX512-float: 16
117 /// - AVX512-double: 8
118 ///
119 /// ### Examples
120 ///
121 /// ~~~{.cpp}
122 /// std::vector<double> numbers(1000);
123 /// for (std::size_t i=0; i<1000; ++i) {
124 /// numbers[i] = rand();
125 /// }
126 ///
127 /// ROOT::Math::KahanSum<double, 4> k;
128 /// k.Add(numbers.begin(), numbers.end());
129 /// // or
130 /// k.Add(numbers);
131 /// ~~~
132 /// ~~~{.cpp}
133 /// double offset = 10.;
134 /// auto result = ROOT::Math::KahanSum<double, 4>::Accumulate(numbers.begin(), numbers.end(), offset);
135 /// ~~~
136 template<typename T = double, unsigned int N = 1>
137 class KahanSum {
138 public:
139 /// Initialise the sum.
140 /// \param[in] initialValue Initialise with this value. Defaults to 0.
141 explicit KahanSum(T initialValue = T{}) {
142 fSum[0] = initialValue;
143 std::fill(std::begin(fSum)+1, std::end(fSum), 0.);
144 std::fill(std::begin(fCarry), std::end(fCarry), 0.);
145 }
146
147 /// Initialise with a sum value and a carry value.
148 /// \param[in] initialSumValue Initialise the sum with this value.
149 /// \param[in] initialCarryValue Initialise the carry with this value.
153 std::fill(std::begin(fSum)+1, std::end(fSum), 0.);
154 std::fill(std::begin(fCarry)+1, std::end(fCarry), 0.);
155 }
156
157 /// Initialise the sum with a pre-existing state.
158 /// \param[in] sumBegin Begin of iterator range with values to initialise the sum with.
159 /// \param[in] sumEnd End of iterator range with values to initialise the sum with.
160 /// \param[in] carryBegin Begin of iterator range with values to initialise the carry with.
161 /// \param[in] carryEnd End of iterator range with values to initialise the carry with.
162 template<class Iterator>
163 KahanSum(Iterator sumBegin, Iterator sumEnd, Iterator carryBegin, Iterator carryEnd) {
164 assert(std::distance(sumBegin, sumEnd) == N);
165 assert(std::distance(carryBegin, carryEnd) == N);
166 std::copy(sumBegin, sumEnd, std::begin(fSum));
167 std::copy(carryBegin, carryEnd, std::begin(fCarry));
168 }
169
170 /// Constructor to create a KahanSum from another KahanSum with a different number of accumulators
171 template <unsigned int M>
173 fSum[0] = other.Sum();
174 fCarry[0] = other.Carry();
175 std::fill(std::begin(fSum)+1, std::end(fSum), 0.);
176 std::fill(std::begin(fCarry)+1, std::end(fCarry), 0.);
177 }
178
179 /// Single-element accumulation. Will not vectorise.
180 void Add(T x) {
181 auto y = x - fCarry[0];
182 auto t = fSum[0] + y;
183 fCarry[0] = (t - fSum[0]) - y;
184 fSum[0] = t;
185 }
186
187
188 /// Accumulate from a range denoted by iterators.
189 ///
190 /// This function will auto-vectorise with random-access iterators.
191 /// \param[in] begin Beginning of a range. Needs to be a random access iterator for automatic
192 /// vectorisation, because a contiguous block of memory needs to be read.
193 /// \param[in] end End of the range.
194 template <class Iterator>
195 void Add(Iterator begin, Iterator end) {
196 static_assert(std::is_floating_point<
197 typename std::remove_reference<decltype(*begin)>::type>::value,
198 "Iterator needs to point to floating-point values.");
199 const std::size_t n = std::distance(begin, end);
200
201 for (std::size_t i=0; i<n; ++i) {
202 AddIndexed(*(begin++), i);
203 }
204 }
205
206
207 /// Fill from a container that supports index access.
208 /// \param[in] inputs Container with index access such as std::vector or array.
209 template<class Container_t>
210 void Add(const Container_t& inputs) {
211 static_assert(std::is_floating_point<typename Container_t::value_type>::value,
212 "Container does not hold floating-point values.");
213 for (std::size_t i=0; i < inputs.size(); ++i) {
214 AddIndexed(inputs[i], i);
215 }
216 }
217
218
219 /// Iterate over a range and return an instance of a KahanSum.
220 ///
221 /// See Add(Iterator,Iterator) for details.
222 /// \param[in] begin Beginning of a range.
223 /// \param[in] end End of the range.
224 /// \param[in] initialValue Optional initial value.
225 template <class Iterator>
226 static KahanSum<T, N> Accumulate(Iterator begin, Iterator end,
227 T initialValue = T{}) {
229 theSum.Add(begin, end);
230
231 return theSum;
232 }
233
234
235 /// Add `input` to the sum.
236 ///
237 /// Particularly helpful when filling from a for loop.
238 /// This function can be inlined and auto-vectorised if
239 /// the index parameter is used to enumerate *consecutive* fills.
240 /// Use Add() or Accumulate() when no index is available.
241 /// \param[in] input Value to accumulate.
242 /// \param[in] index Index of the value. Determines internal accumulator that this
243 /// value is added to. Make sure that consecutive fills have consecutive indices
244 /// to make a loop auto-vectorisable. The actual value of the index does not matter,
245 /// as long as it is consecutive.
246 void AddIndexed(T input, std::size_t index) {
247 const unsigned int i = index % N;
248 const T y = input - fCarry[i];
249 const T t = fSum[i] + y;
250 fCarry[i] = (t - fSum[i]) - y;
251 fSum[i] = t;
252 }
253
254 /// \return Compensated sum.
255 T Sum() const {
256 return std::accumulate(std::begin(fSum), std::end(fSum), 0.);
257 }
258
259 /// \return Compensated sum.
260 T Result() const {
261 return Sum();
262 }
263
264 /// \return The sum used for compensation.
265 T Carry() const {
266 return std::accumulate(std::begin(fCarry), std::end(fCarry), 0.);
267 }
268
269 /// Add `arg` into accumulator. Does not vectorise.
271 Add(arg);
272 return *this;
273 }
274
275 /// Add other KahanSum into accumulator. Does not vectorise.
276 ///
277 /// Based on KahanIncrement from:
278 /// Y. Tian, S. Tatikonda and B. Reinwald, "Scalable and Numerically Stable Descriptive Statistics in SystemML," 2012 IEEE 28th International Conference on Data Engineering, 2012, pp. 1351-1359, doi: 10.1109/ICDE.2012.12.
279 /// Note that while Tian et al. add the carry in the first step, we subtract
280 /// the carry, in accordance with the Add(Indexed) implementation(s) above.
281 /// This is purely an implementation choice that has no impact on performance.
282 ///
283 /// \note Take care when using += (and -=) to add other KahanSums into a zero-initialized
284 /// KahanSum. The operator behaves correctly in this case, but the result may be slightly
285 /// off if you expect 0 + x to yield exactly x (where 0 is the zero-initialized KahanSum
286 /// and x another KahanSum). In particular, x's carry term may get lost. This doesn't
287 /// just happen with zero-initialized KahanSums; see the SubtractWithABitTooSmallCarry
288 /// test case in the testKahan unittest for other examples. This behavior is internally
289 /// consistent: the carry also gets lost if you switch the operands and it also happens with
290 /// other KahanSum operators.
291 template<typename U, unsigned int M>
293 U corrected_arg_sum = other.Sum() - (fCarry[0] + other.Carry());
294 U sum = fSum[0] + corrected_arg_sum;
296 fSum[0] = sum;
297 fCarry[0] = correction;
298 return *this;
299 }
300
301 /// Subtract other KahanSum. Does not vectorise.
302 ///
303 /// Based on KahanIncrement from: Tian et al., 2012 (see operator+= documentation).
304 template<typename U, unsigned int M>
306 U corrected_arg_sum = -other.Sum() - (fCarry[0] - other.Carry());
307 U sum = fSum[0] + corrected_arg_sum;
309 fSum[0] = sum;
310 fCarry[0] = correction;
311 return *this;
312 }
313
315 {
316 return {-this->fSum[0], -this->fCarry[0]};
317 }
318
319 template<typename U, unsigned int M>
320 bool operator ==(KahanSum<U, M> const& other) const {
321 return (this->Sum() == other.Sum()) && (this->Carry() == other.Carry());
322 }
323
324 template<typename U, unsigned int M>
325 bool operator !=(KahanSum<U, M> const& other) const {
326 return !(*this == other);
327 }
328
329 private:
330 T fSum[N];
332 };
333
334 /// Add two non-vectorized KahanSums.
335 template<typename T, unsigned int N, typename U, unsigned int M>
337 KahanSum<T, N> sum(left);
338 sum += right;
339 return sum;
340 }
341
342 /// Subtract two non-vectorized KahanSums.
343 template<typename T, unsigned int N, typename U, unsigned int M>
345 KahanSum<T, N> sum(left);
346 sum -= right;
347 return sum;
348 }
349
350 } // end namespace Math
351
352} // end namespace ROOT
353
354
355#endif /* ROOT_Math_Util */
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define N
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition Util.h:137
T Sum() const
Definition Util.h:255
bool operator==(KahanSum< U, M > const &other) const
Definition Util.h:320
KahanSum< T, N > operator-()
Definition Util.h:314
static KahanSum< T, N > Accumulate(Iterator begin, Iterator end, T initialValue=T{})
Iterate over a range and return an instance of a KahanSum.
Definition Util.h:226
T Result() const
Definition Util.h:260
void Add(Iterator begin, Iterator end)
Accumulate from a range denoted by iterators.
Definition Util.h:195
void Add(const Container_t &inputs)
Fill from a container that supports index access.
Definition Util.h:210
KahanSum(KahanSum< T, M > const &other)
Constructor to create a KahanSum from another KahanSum with a different number of accumulators.
Definition Util.h:172
bool operator!=(KahanSum< U, M > const &other) const
Definition Util.h:325
T Carry() const
Definition Util.h:265
void AddIndexed(T input, std::size_t index)
Add input to the sum.
Definition Util.h:246
KahanSum< T, N > & operator+=(T arg)
Add arg into accumulator. Does not vectorise.
Definition Util.h:270
KahanSum< T, N > & operator+=(const KahanSum< U, M > &other)
Add other KahanSum into accumulator.
Definition Util.h:292
KahanSum(T initialValue=T{})
Initialise the sum.
Definition Util.h:141
KahanSum(T initialSumValue, T initialCarryValue)
Initialise with a sum value and a carry value.
Definition Util.h:150
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition Util.h:180
KahanSum< T, N > & operator-=(KahanSum< U, M > const &other)
Subtract other KahanSum.
Definition Util.h:305
KahanSum(Iterator sumBegin, Iterator sumEnd, Iterator carryBegin, Iterator carryEnd)
Initialise the sum with a pre-existing state.
Definition Util.h:163
std::chrono::steady_clock::time_point fBegin
Definition Util.h:56
const std::string fMessage
Definition Util.h:58
std::function< void(std::string const &) fPrinter)
Definition Util.h:57
TimingScope(std::function< void(std::string const &)> printer, std::string const &message)
Definition Util.cxx:3
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Namespace for new Math classes and functions.
T EvalLog(T x)
safe evaluation of log(x) with a protections against negative or zero argument to the log smooth line...
Definition Util.h:79
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition Util.h:65
DisplacementVector2D< CoordSystem1, U > operator+(DisplacementVector2D< CoordSystem1, U > v1, const DisplacementVector2D< CoordSystem2, U > &v2)
Addition of DisplacementVector2D vectors.
DisplacementVector2D< CoordSystem1, U > operator-(DisplacementVector2D< CoordSystem1, U > v1, DisplacementVector2D< CoordSystem2, U > const &v2)
Difference between two DisplacementVector2D vectors.
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339