Belle II Software light-2406-ragdoll
CalcMeanCov.h
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#pragma once
10
11#include <cmath>
12
13namespace Belle2 {
34 template<int N = 1, class RealType = double>
36 static_assert(N > 0, "Number of parameters, N, must be positive");
37 public:
40
42 typedef RealType value_type;
43
45 void clear()
46 {
47 m_entries = 0;
48 for (int i = 0; i < N; i++)
49 m_mean[i] = 0;
50 for (int i = 0; i < N * (N + 1) / 2; i++)
51 m_covariance[i] = 0;
52 }
53
59 template<class... T> void addWeighted(value_type weight, T... values)
60 {
61 static_assert(sizeof...(values) == N,
62 "Number of arguments must be equal N");
63 m_entries += weight;
64 addValue<0>(weight, values...);
65 }
66
71 template<class... T> void add(T... values)
72 {
73 addWeighted(1.0, values...);
74 }
75
77 void add(const CalcMeanCov<N, RealType>& other)
78 {
79 const value_type n = m_entries + other.m_entries;
80 if (n == 0)
81 return;
82 for (int i = 0; i < N; ++i) {
83 const value_type delta = (other.m_mean[i] - m_mean[i]);
84 //Update covariance matrix
85 for (int j = i; j < N; ++j) {
86 const value_type delta2 = (other.m_mean[j] - m_mean[j]);
87 m_covariance[getIndex(i, j)] += other.m_covariance[getIndex(i, j)] + m_entries * other.m_entries
88 / n * delta * delta2;
89 }
90 //Update mean value
91 m_mean[i] += other.m_entries * delta / n;
92 }
93 m_entries = n;
94 }
95
101 {
102 m_entries += weight;
103 addArrayValues(weight, values);
104 }
105
109 void addArray(value_type* values)
110 {
111 addWeightedArray(1.0, values);
112 }
113
120 value_type getEntries() const { return m_entries; }
122 value_type getMean(int i) const { return m_mean[i]; }
124 value_type getCovariance(int i, int j) const
125 {
126 if (m_entries == 0.0) return 0.0;
127 return m_covariance[getIndex(i, j)] / m_entries;
128 }
130 value_type getCorrelation(int i, int j) const
131 {
132 if (i == j) return 1.0;
133 if (m_entries == 0.0) return 0.0;
134 return getCovariance(i, j) / (getStddev(i) * getStddev(j));
135 }
137 value_type getVariance(int i) const { return getCovariance(i, i); }
139 value_type getStddev(int i) const { return std::sqrt(getVariance(i)); }
141 value_type getSum(int i) const { return getMean(i) * m_entries; }
142
152 template <int i = 0> value_type getMean() const
153 {
154 static_assert(i >= 0 && i < N, "index i out of range");
155 return m_mean[i];
156 }
158 template <int i, int j> value_type getCovariance() const
159 {
160 static_assert(i >= 0 && i < N, "index i out of range");
161 static_assert(j >= 0 && j < N, "index j out of range");
162 if (m_entries == 0.0) return 0.0;
163 return m_covariance[getIndex(i, j)] / m_entries;
164 }
166 template <int i, int j> value_type getCorrelation() const
167 {
168 if (i == j) return 1.0;
169 if (m_entries == 0.0) return 0.0;
170 return getCovariance<i, j>() / (getStddev<i>() * getStddev<j>());
171 }
173 template <int i = 0> value_type getVariance() const
174 {
175 return getCovariance<i, i>();
176 }
178 template <int i = 0> value_type getStddev() const
179 {
180 return std::sqrt(getVariance<i>());
181 }
183 template <int i = 0> value_type getSum() const
184 {
185 return getMean<i>() * m_entries;
186 }
187
190 private:
198 template<int i, class... T>
199 void addValue(value_type weight, value_type x, T... values)
200 {
201 const value_type delta = (x - m_mean[i]);
202 //Update covariance elements
203 updateCov<i, i>(weight, delta, x, values...);
204 //Update mean value
205 m_mean[i] += weight * delta / m_entries;
206 //recursively call again for remaining values
207 addValue < i + 1 > (weight, values...);
208 }
210 template<int i> void addValue(value_type) {}
211
221 template<int i, int j, class... T>
223 T... values)
224 {
225 const value_type delta2 = (x - m_mean[j]);
226 m_covariance[getIndex(i, j)] += (m_entries - weight) * weight
227 / m_entries * delta * delta2;
228 updateCov < i, j + 1 > (weight, delta, values...);
229 }
230
232 template<int i, int j> void updateCov(value_type, value_type) {}
233
241 {
242 for (int i = 0; i < N; ++i) {
243 const value_type delta = (x[i] - m_mean[i]);
244 //Update covariance matrix
245 for (int j = i; j < N; ++j) {
246 const value_type delta2 = (x[j] - m_mean[j]);
247 m_covariance[getIndex(i, j)] += (m_entries - weight) * weight
248 / m_entries * delta * delta2;
249 }
250 //Update mean value
251 m_mean[i] += weight * delta / m_entries;
252 }
253 }
254
260 constexpr int getIndex(unsigned int i, unsigned int j) const
261 {
262 //swap indices if i >= j
263 return (i < j) ? ((j + 1) * j / 2 + i) : ((i + 1) * i / 2 + j);
264 }
265
272 value_type m_covariance[N * (N + 1) / 2];
273 };
274
276} //Belle2 namespace
Class to calculate mean and and covariance between a number of parameters on running data without sto...
Definition: CalcMeanCov.h:35
value_type m_entries
Store the sum of weights.
Definition: CalcMeanCov.h:267
void addWeighted(value_type weight, T... values)
Update mean and covariance by adding a new, weighted entry.
Definition: CalcMeanCov.h:59
value_type getVariance(int i) const
Return the variance for paramter i.
Definition: CalcMeanCov.h:137
void add(const CalcMeanCov< N, RealType > &other)
Merge the data set in 'other' into this one.
Definition: CalcMeanCov.h:77
void addValue(value_type)
Break recursion of addValue when no parameters are left.
Definition: CalcMeanCov.h:210
void updateCov(value_type weight, value_type delta, value_type x, T... values)
Update covariance between parameters i and j.
Definition: CalcMeanCov.h:222
value_type getCorrelation(int i, int j) const
Return the correlation coefficient between parameters i and j.
Definition: CalcMeanCov.h:130
value_type getStddev() const
Return the standard deviation for parameter i.
Definition: CalcMeanCov.h:178
value_type getStddev(int i) const
Return the standard deviation for parameter i.
Definition: CalcMeanCov.h:139
value_type getMean(int i) const
Return the mean for parameter i.
Definition: CalcMeanCov.h:122
value_type getMean() const
Return the mean for parameter i.
Definition: CalcMeanCov.h:152
constexpr int getIndex(unsigned int i, unsigned int j) const
Access element in triangular matrix including diagonal elements.
Definition: CalcMeanCov.h:260
RealType value_type
type of float variable to use for calculations and storage
Definition: CalcMeanCov.h:42
value_type getVariance() const
Return the variance for paramter i.
Definition: CalcMeanCov.h:173
value_type m_mean[N]
Store the mean values for all parameters.
Definition: CalcMeanCov.h:269
value_type getSum() const
Return the weighted sum values for parameter i.
Definition: CalcMeanCov.h:183
value_type getCovariance(int i, int j) const
Return the covariance between parameters i and j.
Definition: CalcMeanCov.h:124
value_type getCovariance() const
Return the covariance between parameters i and j.
Definition: CalcMeanCov.h:158
value_type m_covariance[N *(N+1)/2]
Store the triangular covariance matrix for all parameters in continous memory.
Definition: CalcMeanCov.h:272
void addArrayValues(value_type weight, value_type *x)
Add a new set of values and update mean and covariance.
Definition: CalcMeanCov.h:240
void addArray(value_type *values)
Update mean and covariance by adding a new entry.
Definition: CalcMeanCov.h:109
value_type getCorrelation() const
Return the correlation coefficient between parameters i and j.
Definition: CalcMeanCov.h:166
void addWeightedArray(value_type weight, value_type *values)
Update mean and covarianced by adding a new weighted entry.
Definition: CalcMeanCov.h:100
void updateCov(value_type, value_type)
Break recursion of updateCov once all parameters are consumed.
Definition: CalcMeanCov.h:232
value_type getSum(int i) const
Return the weighted sum values for parameter i.
Definition: CalcMeanCov.h:141
CalcMeanCov()
default constructor.
Definition: CalcMeanCov.h:39
void add(T... values)
Update mean and covariance by adding a new entry.
Definition: CalcMeanCov.h:71
void clear()
Clear all values.
Definition: CalcMeanCov.h:45
void addValue(value_type weight, value_type x, T... values)
Add a single value for parameter i and update mean and covariance.
Definition: CalcMeanCov.h:199
value_type getEntries() const
Return the number of entries.
Definition: CalcMeanCov.h:120
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24