File indexing completed on 2025-04-19 09:13:40
0001
0002
0003
0004
0005
0006 #ifndef YODA_Dbn_h
0007 #define YODA_Dbn_h
0008
0009 #include "YODA/Exceptions.h"
0010 #include "YODA/Utils/DbnUtils.h"
0011 #include "YODA/Utils/MathUtils.h"
0012 #include "YODA/Utils/MetaUtils.h"
0013 #include <cmath>
0014 #include <string>
0015 #include <array>
0016
0017 namespace YODA {
0018
0019
0020
0021
0022 namespace {
0023
0024
0025 template <size_t ... Is>
0026 constexpr bool has_match(const size_t i, std::integer_sequence<size_t, Is...>) {
0027 return ((i == Is) || ...);
0028 }
0029
0030
0031 template<typename INT>
0032 constexpr std::integer_sequence<INT>
0033 concat_sequences(std::integer_sequence<INT>){
0034 return {};
0035 }
0036
0037
0038 template<typename INT, INT... s, INT... t>
0039 constexpr std::integer_sequence<INT,s...,t...>
0040 concat_sequences(std::integer_sequence<INT, s...>, std::integer_sequence<INT, t...>){
0041 return {};
0042 }
0043
0044
0045 template<typename INT, INT... s, INT... t, class... R>
0046 constexpr auto
0047 concat_sequences(std::integer_sequence<INT, s...>, std::integer_sequence<INT, t...>, R...) {
0048 return concat_sequences(std::integer_sequence<INT,s...,t...>{}, R{}...);
0049 }
0050
0051
0052
0053 template<class INT, INT a, class CHECK>
0054 constexpr auto AxisFilterSingle(std::integer_sequence<INT, a>, CHECK pass) {
0055 if constexpr (pass(a))
0056 return std::integer_sequence<INT, a>{};
0057 else
0058 return std::integer_sequence<INT>{};
0059 }
0060
0061
0062
0063
0064 template<class INT, INT... b, class CHECK>
0065 constexpr auto AxisFilter(std::integer_sequence<INT, b...>, [[maybe_unused]] CHECK pass) {
0066 if constexpr (sizeof...(b) > 0)
0067 return concat_sequences(AxisFilterSingle(std::integer_sequence<INT, b>{}, pass)...);
0068 else
0069 return std::integer_sequence<INT>{};
0070 }
0071
0072 }
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085 template <size_t N>
0086 class DbnBase {
0087 public:
0088
0089 using BaseT = DbnBase<N>;
0090 using FillDim = std::integral_constant<size_t, N>;
0091 using DataSize = std::integral_constant<size_t, 1 + 2*(N+1) + (N*(N-1)/2)>;
0092
0093
0094
0095
0096
0097 DbnBase() { reset(); }
0098
0099
0100
0101
0102
0103
0104 template<size_t dim = N, typename = std::enable_if_t<(dim < 2)>>
0105 DbnBase(const double numEntries, const std::array<double,N+1>& sumW, const std::array<double,N+1>& sumW2)
0106 : _numEntries(numEntries), _sumW(sumW), _sumW2(sumW2) { }
0107
0108
0109
0110
0111
0112 template<size_t dim = N, typename = std::enable_if_t<(dim >= 2)>>
0113 DbnBase(const double numEntries, const std::array<double,N+1>& sumW, const std::array<double,N+1>& sumW2,
0114 const std::array<double,N*(N-1)/2>& sumWcross)
0115 : _numEntries(numEntries), _sumW(sumW), _sumW2(sumW2), _sumWcross(sumWcross) { }
0116
0117
0118
0119
0120 DbnBase(const DbnBase&) = default;
0121
0122
0123
0124
0125 DbnBase(DbnBase&&) = default;
0126
0127
0128
0129
0130 DbnBase& operator = (const DbnBase&) = default;
0131
0132
0133
0134
0135 DbnBase& operator = (DbnBase&&) = default;
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152 void fill(const std::array<double, N> vals, const double weight = 1.0, const double fraction = 1.0) {
0153 fill(vals, weight, fraction, std::make_index_sequence<N>{});
0154 }
0155
0156
0157
0158
0159
0160
0161 template <typename... Args>
0162 void fill(Args&&... args) {
0163
0164 static_assert(sizeof...(args) <= N + 2,
0165 "Value's dimension doesn't match distribution's dimension.");
0166 std::array<double, sizeof...(args)> vals = {{args...}};
0167 double weight = 1.0;
0168 double fraction = 1.0;
0169 if (vals.size() > N) {
0170 weight = vals[N];
0171 if (vals.size() > N + 1) fraction = vals[N + 1];
0172 }
0173 const double sf = fraction * weight;
0174 _numEntries += fraction;
0175 _sumW.at(0) += sf;
0176 _sumW2.at(0) += fraction*sqr(weight);
0177 for (unsigned int i = 0; i < N; ++i) {
0178 _sumW.at(i + 1) += sf * vals.at(i);
0179 _sumW2.at(i + 1) += sf * sqr(vals.at(i));
0180 }
0181
0182 if constexpr(N >= 2) {
0183 size_t idx = 0;
0184 for (size_t i = 0; i < (N-1); ++i) {
0185 for (size_t j = i+1; j < N; ++j) {
0186 _sumWcross.at(idx++) += sf * vals.at(i) * vals.at(j);
0187 }
0188 }
0189 }
0190
0191 }
0192
0193
0194
0195
0196
0197 void set(const double numEntries, const std::array<double,N+1>& sumW,
0198 const std::array<double,N+1>& sumW2,
0199 const std::array<double,N*(N-1)/2>& sumWcross) {
0200 _numEntries = numEntries;
0201 _sumW = sumW;
0202 _sumW2 = sumW2;
0203 _sumWcross = sumWcross;
0204 }
0205
0206
0207
0208
0209 void set(const double numEntries, const std::vector<double>& sumW,
0210 const std::vector<double>& sumW2,
0211 const std::vector<double>& sumWcross = {}) {
0212 if (!(sumW.size() <= (N + 1) || sumW2.size() <= (N + 1) || sumWcross.size() <= (N*(N-1)/2)))
0213 throw UserError("Value's dimension doesn't match distribution's dimension.");
0214 _numEntries = numEntries;
0215 std::copy_n(std::make_move_iterator(sumW.begin()), sumW.size(), _sumW.begin());
0216 std::copy_n(std::make_move_iterator(sumW2.begin()), sumW2.size(), _sumW2.begin());
0217 std::copy_n(std::make_move_iterator(sumWcross.begin()), sumWcross.size(), _sumWcross.begin());
0218 }
0219
0220
0221
0222 void set(const DbnBase<N>& other) {
0223 if (this != &other) *this = other;
0224 }
0225
0226
0227 void set(DbnBase<N>&& other) {
0228 if (this != &other) *this = std::move(other);
0229 }
0230
0231
0232 void reset() {
0233 _numEntries= 0;
0234 _sumW.fill(0);
0235 _sumW2.fill(0);
0236 _sumWcross.fill(0);
0237 }
0238
0239
0240
0241 void scaleW(const double scalefactor) {
0242 _sumW.at(0) *= scalefactor;
0243 _sumW2.at(0) *= sqr(scalefactor);
0244 for (size_t i = 0; i< N; ++i) {
0245
0246 _sumW.at(i+1) *= scalefactor;
0247 _sumW2.at(i+1) *= scalefactor;
0248 }
0249 for (size_t i = 0; i < _sumWcross.size(); ++i) {
0250 _sumWcross.at(i) *= scalefactor;
0251 }
0252 }
0253
0254 void scale(const size_t dim, const double factor) {
0255 if (dim >= N)
0256 throw RangeError("Dimension index should be less than "+std::to_string(N));
0257 _sumW.at(dim+1) *= factor;
0258 _sumW2.at(dim+1) *= sqr(factor);
0259 size_t idx = 0;
0260 for (size_t i = 0; i < (N-1); ++i) {
0261 for (size_t j = i+1; j < N; ++j) {
0262 if (i == dim || j == dim) {
0263 _sumWcross.at(idx++) *= factor;
0264 }
0265 }
0266 }
0267 }
0268
0269
0270
0271
0272 public:
0273
0274
0275
0276
0277
0278
0279 double errW() const;
0280
0281
0282 double relErrW() const;
0283
0284
0285 double mean(const size_t i) const;
0286
0287
0288 double variance(const size_t i) const;
0289
0290
0291 double stdDev(const size_t i) const { return std::sqrt(variance(i)); }
0292
0293
0294 double stdErr(const size_t i) const;
0295
0296
0297 double relStdErr(const size_t i) const;
0298
0299
0300 double RMS(const size_t i) const;
0301
0302
0303
0304
0305
0306
0307
0308
0309 double numEntries() const {
0310 return _numEntries;
0311 }
0312
0313
0314 double effNumEntries() const {
0315 if (_sumW2.at(0) == 0) return 0;
0316 return _sumW.at(0)*_sumW.at(0) / _sumW2.at(0);
0317 }
0318
0319
0320 double sumW() const {
0321 return _sumW.at(0);
0322 }
0323
0324
0325 double sumW2() const {
0326 return _sumW2.at(0);
0327 }
0328
0329
0330 double sumW(const size_t i) const {
0331 return _sumW.at(i);
0332 }
0333
0334
0335 double sumW2(const size_t i) const {
0336 return _sumW2.at(i);
0337 }
0338
0339
0340 template<size_t dim = N, typename = std::enable_if_t<(dim >= 2)>>
0341 double crossTerm(const size_t A1, const size_t A2) const {
0342 if (A1 >= N || A2 >= N) throw RangeError("Invalid axis int, must be in range 0..dim-1");
0343 if (A1 >= A2) throw RangeError("Indices need to be different for cross term");
0344
0345 size_t idx = 0;
0346 for (size_t i = 0; i < (N-1); ++i) {
0347 for (size_t j = i+1; j < N; ++j) {
0348 if (i == A1 && j == A2) break;
0349 ++idx;
0350 }
0351 if (i == A1) break;
0352 }
0353
0354 return _sumWcross.at(idx);
0355 }
0356
0357 size_t dim() const { return N; }
0358
0359
0360
0361
0362
0363
0364
0365 std::string toString() const {
0366 std::string res ="";
0367 res += ("numEntries="+ std::to_string(_numEntries)) ;
0368 if (sumW()) {
0369 res += (", Mean="+ std::to_string(mean(0))) ;
0370 res += (", RMS="+ std::to_string(RMS(0))) ;
0371 }
0372 return res;
0373 }
0374
0375
0376
0377
0378
0379
0380
0381
0382 template<typename... size_t>
0383 DbnBase<sizeof...(size_t)> reduceTo(size_t... axes) const {
0384 if constexpr (N == sizeof...(axes)) { return *this; }
0385 else {
0386
0387 std::array<double, sizeof...(axes)+1> newSumW = { _sumW[0], _sumW[axes+1]... };
0388 std::array<double, sizeof...(axes)+1> newSumW2 = { _sumW2[0], _sumW2[axes+1]... };
0389 if constexpr (sizeof...(axes) < 2) {
0390 return DbnBase<sizeof...(axes)>(_numEntries, newSumW, newSumW2);
0391 }
0392 else {
0393 constexpr auto newDim = sizeof...(axes);
0394 std::array<double, newDim*(newDim-1)/2> newSumWcross;
0395 unsigned int idx = 0;
0396 for (unsigned int i : {axes...}) {
0397 for (unsigned int j : {axes...}) {
0398 if (i < j) { newSumWcross.at(idx) = crossTerm(i, j); ++idx; }
0399 }
0400 }
0401 return DbnBase<sizeof...(axes)>(_numEntries, newSumW, newSumW2, newSumWcross);
0402 }
0403 }
0404 }
0405
0406
0407 template<size_t... Is>
0408 DbnBase<sizeof...(Is)> reduceTo(std::index_sequence<Is...>) const {
0409 return reduceTo(Is...);
0410 }
0411
0412
0413
0414 template<size_t... axes>
0415 auto reduce() const {
0416 constexpr auto new_axes = AxisFilter(std::make_index_sequence<N>{},
0417 [](size_t i) { return !has_match(i, std::integer_sequence<size_t, axes...>{}); }
0418 );
0419 return reduceTo(new_axes);
0420 }
0421
0422
0423
0424
0425
0426
0427 size_t _lengthContent() const noexcept {
0428 return DataSize::value;
0429 }
0430
0431 std::vector<double> _serializeContent() const noexcept {
0432 std::vector<double> rtn;
0433 rtn.reserve(DataSize::value);
0434 rtn.insert(rtn.end(), _sumW.begin(), _sumW.end());
0435 rtn.insert(rtn.end(), _sumW2.begin(), _sumW2.end());
0436 rtn.insert(rtn.end(), _sumWcross.begin(), _sumWcross.end());
0437 rtn.push_back(_numEntries);
0438 return rtn;
0439 }
0440
0441 void _deserializeContent(const std::vector<double>& data) {
0442
0443 if (data.size() != DataSize::value)
0444 throw UserError("Length of serialized data should be "+std::to_string(DataSize::value)+"!");
0445
0446 auto itr = data.cbegin();
0447 std::copy_n(itr, _sumW.size(), _sumW.begin());
0448 std::copy_n(itr+N+1, _sumW2.size(), _sumW2.begin());
0449 std::copy_n(itr+2*(N+1), _sumWcross.size(), _sumWcross.begin());
0450 _numEntries = *(itr + 2*(N+1) + (N*(N-1)/2));
0451
0452 }
0453
0454
0455
0456
0457
0458
0459
0460 DbnBase& operator += (const DbnBase& d);
0461
0462
0463 DbnBase& operator += (DbnBase&& d);
0464
0465
0466 DbnBase& operator -= (const DbnBase& d);
0467
0468
0469 DbnBase& operator -= (DbnBase&& d);
0470
0471
0472
0473 private:
0474
0475
0476
0477
0478 double _numEntries;
0479
0480
0481 std::array<double, N+1> _sumW;
0482
0483
0484 std::array<double, N+1> _sumW2;
0485
0486
0487 std::array<double, N*(N-1)/2> _sumWcross;
0488
0489
0490
0491
0492
0493
0494
0495 template <size_t... I>
0496 void fill(const std::array<double, N> vals, const double weight, const double fraction,
0497 std::index_sequence<I...>) {
0498 fill(vals[I]..., weight, fraction);
0499 }
0500
0501
0502 };
0503
0504
0505 template <size_t N>
0506 double DbnBase<N>::errW() const {
0507 return sqrt(sumW2(0));
0508 }
0509
0510 template <size_t N>
0511 double DbnBase<N>::relErrW() const {
0512 if (effNumEntries() == 0) {
0513 return std::numeric_limits<double>::quiet_NaN();
0514 }
0515 return errW()/sumW(0);
0516 }
0517
0518 template <size_t N>
0519 double DbnBase<N>::mean(const size_t i) const {
0520 return YODA::mean(sumW(i), sumW(0));
0521 }
0522
0523
0524 template <size_t N>
0525 double DbnBase<N>::variance(const size_t i) const {
0526 return YODA::variance(sumW(i), sumW(0), sumW2(i), sumW2(0));
0527 }
0528
0529
0530 template <size_t N>
0531 double DbnBase<N>::stdErr(const size_t i) const {
0532 return YODA::stdErr(sumW(i), sumW(0), sumW2(i), sumW2(0));
0533 }
0534
0535
0536 template <size_t N>
0537 double DbnBase<N>::relStdErr(const size_t i) const {
0538 if (effNumEntries() == 0) {
0539 return std::numeric_limits<double>::quiet_NaN();
0540 }
0541 return stdErr(i) / mean(i);
0542 }
0543
0544
0545 template <size_t N>
0546 double DbnBase<N>::RMS(const size_t i) const {
0547 return YODA::RMS(sumW2(i), sumW(0), sumW2());
0548 }
0549
0550 template <size_t N>
0551 DbnBase<N>& DbnBase<N>::operator += (const DbnBase<N>& d) {
0552 _numEntries += d._numEntries;
0553 for (size_t i = 0; i <= N; ++i) {
0554 _sumW.at(i) += d._sumW.at(i);
0555 _sumW2.at(i) += d._sumW2.at(i);
0556 }
0557 for (size_t i = 0; i < _sumWcross.size(); ++i) {
0558 _sumWcross.at(i) += d._sumWcross.at(i);
0559 }
0560 return *this;
0561 }
0562
0563
0564 template <size_t N>
0565 DbnBase<N>& DbnBase<N>::operator += (DbnBase<N>&& d) {
0566 _numEntries += std::move(d._numEntries);
0567 for (size_t i = 0; i <= N; ++i) {
0568 _sumW.at(i) += std::move(d._sumW.at(i));
0569 _sumW2.at(i) += std::move(d._sumW2.at(i));
0570 }
0571 for (size_t i = 0; i < _sumWcross.size(); ++i) {
0572 _sumWcross.at(i) += std::move(d._sumWcross.at(i));
0573 }
0574 return *this;
0575 }
0576
0577
0578 template <size_t N>
0579 DbnBase<N>& DbnBase<N>::operator -= (const DbnBase<N>& d) {
0580 _numEntries -= d._numEntries;
0581 for (unsigned int i =0; i<= N ; ++i) {
0582 _sumW.at(i) -= d._sumW.at(i);
0583 _sumW2.at(i) += d._sumW2.at(i);
0584 }
0585 for (size_t i = 0; i < _sumWcross.size(); ++i) {
0586 _sumWcross.at(i) -= d._sumWcross.at(i);
0587 }
0588 return *this;
0589 }
0590
0591
0592 template <size_t N>
0593 DbnBase<N>& DbnBase<N>::operator -= (DbnBase<N>&& d) {
0594 _numEntries -= std::move(d._numEntries);
0595 for (unsigned int i =0; i<= N ; ++i) {
0596 _sumW.at(i) -= std::move(d._sumW.at(i));
0597 _sumW2.at(i) += std::move(d._sumW2.at(i));
0598 }
0599 for (size_t i = 0; i < _sumWcross.size(); ++i) {
0600 _sumWcross.at(i) -= std::move(d._sumWcross.at(i));
0601 }
0602 return *this;
0603 }
0604
0605
0606
0607 template <size_t N>
0608 inline DbnBase<N> operator + (DbnBase<N> first, const DbnBase<N>& second) {
0609 first += second;
0610 return first;
0611 }
0612
0613
0614 template <size_t N>
0615 inline DbnBase<N> operator + (DbnBase<N> first, DbnBase<N>&& second) {
0616 first += std::move(second);
0617 return first;
0618 }
0619
0620
0621 template <size_t N>
0622 inline DbnBase<N> operator - (DbnBase<N> first, const DbnBase<N>& second) {
0623 first -= second;
0624 return first;
0625 }
0626
0627
0628 template <size_t N>
0629 inline DbnBase<N> operator - (DbnBase<N> first, DbnBase<N>&& second) {
0630 first -= std::move(second);
0631 return first;
0632 }
0633
0634
0635
0636 template<size_t N>
0637 class Dbn : public DbnBase<N> {
0638 public:
0639 using BaseT = DbnBase<N>;
0640 using BaseT::BaseT;
0641 using BaseT::operator=;
0642
0643 };
0644
0645
0646 template<>
0647 class Dbn<0> : public DbnBase<0> {
0648 public:
0649
0650 using BaseT = DbnBase<0>;
0651 using BaseT::BaseT;
0652 using BaseT::operator=;
0653
0654 Dbn() : BaseT() {}
0655
0656 Dbn(const double numEntries, const double sumW, const double sumW2)
0657 : BaseT(numEntries, {sumW} , {sumW2}) { }
0658
0659 Dbn(const BaseT& other) : BaseT(other) {}
0660
0661 Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0662
0663 void fill(const std::array<double, 0>& vals, const double weight = 1.0, const double fraction = 1.0) {
0664 BaseT::fill(vals, weight, fraction);
0665 }
0666
0667 void fill(const double weight=1.0, const double fraction=1.0) { BaseT::fill({}, weight, fraction); }
0668
0669 void set(const double numEntries, const double sumW, const double sumW2) {
0670 BaseT::set(numEntries, {sumW}, {sumW2});
0671 }
0672
0673 };
0674
0675
0676
0677 template<>
0678 class Dbn<1> : public DbnBase<1>,
0679 public XDbnMixin<Dbn<1>> {
0680 public:
0681
0682 using BaseT = DbnBase<1>;
0683 using BaseT::BaseT;
0684 using BaseT::operator=;
0685
0686 Dbn() : BaseT() {}
0687
0688 Dbn(const double numEntries, const double sumW, const double sumW2, const double sumWX, const double sumWX2)
0689 : BaseT(numEntries, {sumW, sumWX} , {sumW2, sumWX2}) { }
0690
0691 Dbn(const BaseT& other) : BaseT(other) {}
0692
0693 Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0694
0695 void fill(const std::array<double, 1>& vals, const double weight = 1.0, const double fraction = 1.0) {
0696 BaseT::fill(vals, weight, fraction);
0697 }
0698
0699 void fill(const double val, const double weight=1.0, const double fraction=1.0) {
0700 BaseT::fill({val}, weight, fraction);
0701 }
0702
0703 };
0704
0705
0706
0707 template<>
0708 class Dbn<2> : public DbnBase<2>,
0709 public XDbnMixin<Dbn<2>>,
0710 public YDbnMixin<Dbn<2>> {
0711 public:
0712
0713 using BaseT = DbnBase<2>;
0714 using BaseT::BaseT;
0715 using BaseT::operator=;
0716
0717 Dbn() : BaseT() {}
0718
0719 Dbn(const double numEntries, const double sumW, const double sumW2, const double sumWX, const double sumWX2,
0720 const double sumWY, const double sumWY2, const double sumWXY)
0721 : BaseT(numEntries, {sumW, sumWX, sumWY} , {sumW2, sumWX2, sumWY2}, {sumWXY}) { }
0722
0723 Dbn(const BaseT& other) : BaseT(other) {}
0724
0725 Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0726
0727 void fill(const std::array<double, 2>& vals, const double weight = 1.0, const double fraction = 1.0) {
0728 BaseT::fill(vals, weight, fraction);
0729 }
0730
0731 void fill(const double valX, const double valY, const double weight=1.0, const double fraction=1.0) {
0732 BaseT::fill({valX, valY}, weight, fraction);
0733 }
0734
0735
0736 };
0737
0738
0739
0740 template<>
0741 class Dbn<3> : public DbnBase<3>,
0742 public XDbnMixin<Dbn<3>>,
0743 public YDbnMixin<Dbn<3>>,
0744 public ZDbnMixin<Dbn<3>> {
0745 public:
0746
0747 using BaseT = DbnBase<3>;
0748 using BaseT::BaseT;
0749 using BaseT::operator=;
0750
0751 Dbn() : BaseT() {}
0752
0753 Dbn(const double numEntries, const double sumW, const double sumW2,
0754 const double sumWX, const double sumWX2,
0755 const double sumWY, const double sumWY2,
0756 const double sumWZ, const double sumWZ2,
0757 const double sumWXY, const double sumWXZ, const double sumWYZ)
0758 : BaseT(numEntries, {sumW, sumWX, sumWY, sumWZ} , {sumW2, sumWX2, sumWY2, sumWZ2}, {sumWXY, sumWXZ, sumWYZ}) { }
0759
0760 Dbn(const BaseT& other) : BaseT(other) {}
0761
0762 Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0763
0764 void fill(const std::array<double, 3>& vals, const double weight = 1.0, const double fraction = 1.0) {
0765 BaseT::fill(vals, weight, fraction);
0766 }
0767
0768 void fill(const double valX, const double valY, const double valZ, const double weight=1.0, const double fraction=1.0) {
0769 BaseT::fill({valX, valY, valZ}, weight, fraction);
0770 }
0771
0772 };
0773
0774
0775 using Dbn0D = Dbn<0>;
0776 using Dbn1D = Dbn<1>;
0777 using Dbn2D = Dbn<2>;
0778 using Dbn3D = Dbn<3>;
0779 }
0780
0781 #endif