ROOT   Reference Guide
df022_useKahan.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_dataframe
3/// \notebook
4/// Implement a custom action that evaluates a Kahan sum.
5///
6/// This tutorial shows how to implement a Kahan summation custom action.
7///
8/// \macro_code
9/// \macro_output
10///
11/// \date July 2018
12/// \authors Enrico Guiraud, Danilo Piparo (CERN), Massimo Tumolo (Politecnico di Torino)
13
14template <typename T>
15class KahanSum final : public ROOT::Detail::RDF::RActionImpl<class KahanSum<T>> {
16public:
17 /// This type is a requirement for every helper.
18 using Result_t = T;
19
20private:
21 std::vector<T> fPartialSums;
22 std::vector<T> fCompensations;
23 int fNSlots;
24
25 std::shared_ptr<T> fResultSum;
26
27 void KahanAlgorithm(const T &x, T &sum, T &compensation){
28 T y = x - compensation;
29 T t = sum + y;
30 compensation = (t - sum) - y;
31 sum = t;
32 }
33
34public:
35 KahanSum(KahanSum &&) = default;
36 KahanSum(const KahanSum &) = delete;
37
38 KahanSum(const std::shared_ptr<T> &r) : fResultSum(r)
39 {
40 static_assert(std::is_floating_point<T>::value, "Kahan sum makes sense only on floating point numbers");
41
43 fPartialSums.resize(fNSlots, 0.);
44 fCompensations.resize(fNSlots, 0.);
45 }
46
47 std::shared_ptr<Result_t> GetResultPtr() const { return fResultSum; }
48
49 void Initialize() {}
50 void InitTask(TTreeReader *, unsigned int) {}
51
52 void Exec(unsigned int slot, T x)
53 {
54 KahanAlgorithm(x, fPartialSums[slot], fCompensations[slot]);
55 }
56
57 template <typename V=T, std::enable_if_t<ROOT::Internal::RDF::IsDataContainer<V>::value, int> = 0>
58 void Exec(unsigned int slot, const T &vs)
59 {
60 for (auto &&v : vs) {
61 Exec(slot, v);
62 }
63 }
64
65 void Finalize()
66 {
67 T sum(0) ;
68 T compensation(0);
69 for (int i = 0; i < fNSlots; ++i) {
70 KahanAlgorithm(fPartialSums[i], sum, compensation);
71 }
72 *fResultSum = sum;
73 }
74
75 std::string GetActionName(){
76 return "THnHelper";
77 }
78
79};
80
81void df022_useKahan()
82{
83 // We enable implicit parallelism
85
87 auto dd = d.Define("x", "(rdfentry_ %2 == 0) ? 0.00000001 : 100000000.");
88
89 auto ptr = std::make_shared<double>();
90 KahanSum<double> helper(ptr);
91
92 auto kahanResult = dd.Book<double>(std::move(helper), {"x"});
93 auto plainResult = dd.Sum<double>({"x"});
94
95 std::cout << std::setprecision(24) << "Kahan: " << *kahanResult << " Classical: " << *plainResult << std::endl;
96 // Outputs: Kahan: 1000000000.00000011920929 Classical: 1000000000
97}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
ROOT's RDataFrame offers a high level interface for analyses of data stored in TTree,...
Definition: RDataFrame.hxx:42
A simple, robust and fast interface to read values from ROOT columnar datasets such as TTree,...
Definition: TTreeReader.h:44
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
CPYCPPYY_EXTERN bool Exec(const std::string &cmd)
Definition: API.cxx:333
double T(double x)
Definition: ChebyshevPol.h:34
void EnableImplicitMT(UInt_t numthreads=0)
Enable ROOT's implicit multi-threading for all objects and methods that provide an internal paralleli...
Definition: TROOT.cxx:524
Bool_t IsImplicitMTEnabled()
Returns true if the implicit multi-threading in ROOT is enabled.
Definition: TROOT.cxx:555
UInt_t GetThreadPoolSize()
Returns the size of ROOT's thread pool.
Definition: TROOT.cxx:562
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.cxx:176
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345