Dbn.h

Go to the documentation of this file.

// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2024 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_Dbn_h
#define YODA_Dbn_h
 
#include "YODA/Exceptions.h"
#include "YODA/Utils/DbnUtils.h"
#include "YODA/Utils/MathUtils.h"
#include "YODA/Utils/MetaUtils.h"
#include <cmath>
#include <string>
#include <array>
 
namespace YODA {
 
  namespace {
 
    template <size_t ... Is>
    constexpr bool has_match(const size_t i, std::integer_sequence<size_t, Is...>) {
      return ((i == Is) || ...);
    }
 
    template<typename INT>
    constexpr std::integer_sequence<INT>
    concat_sequences(std::integer_sequence<INT>){
        return {};
    }
 
    template<typename INT, INT... s, INT... t>
    constexpr std::integer_sequence<INT,s...,t...>
    concat_sequences(std::integer_sequence<INT, s...>, std::integer_sequence<INT, t...>){
        return {};
    }
 
    template<typename INT, INT... s, INT... t, class... R>
    constexpr auto
    concat_sequences(std::integer_sequence<INT, s...>, std::integer_sequence<INT, t...>, R...) {
        return concat_sequences(std::integer_sequence<INT,s...,t...>{}, R{}...);
    }
 
    template<class INT, INT a, class CHECK>
    constexpr auto AxisFilterSingle(std::integer_sequence<INT, a>, CHECK pass) {
      if constexpr (pass(a))
        return std::integer_sequence<INT, a>{};
      else
        return std::integer_sequence<INT>{};
    }
 
    template<class INT, INT... b, class CHECK>
    constexpr auto AxisFilter(std::integer_sequence<INT, b...>, [[maybe_unused]] CHECK pass) {
      if constexpr (sizeof...(b) > 0) // non empty sequence
        return concat_sequences(AxisFilterSingle(std::integer_sequence<INT, b>{}, pass)...);
      else // empty sequence case
        return std::integer_sequence<INT>{};
    }
 
  }
 
  template <size_t N>
  class DbnBase {
  public:
 
    using BaseT = DbnBase<N>;
    using FillDim = std::integral_constant<size_t, N>;
    using DataSize = std::integral_constant<size_t, 1 + 2*(N+1) + (N*(N-1)/2)>;
 
 
 
    DbnBase() { reset(); }
 
 
    template<size_t dim = N, typename = std::enable_if_t<(dim < 2)>>
    DbnBase(const double numEntries, const std::array<double,N+1>& sumW, const std::array<double,N+1>& sumW2)
      : _numEntries(numEntries), _sumW(sumW), _sumW2(sumW2) { }
 
    template<size_t dim = N, typename = std::enable_if_t<(dim >= 2)>>
    DbnBase(const double numEntries, const std::array<double,N+1>& sumW, const std::array<double,N+1>& sumW2,
                                                                         const std::array<double,N*(N-1)/2>& sumWcross)
      : _numEntries(numEntries), _sumW(sumW), _sumW2(sumW2), _sumWcross(sumWcross) { }
 
    DbnBase(const DbnBase&) = default;
 
    DbnBase(DbnBase&&) = default;
 
    DbnBase& operator = (const DbnBase&) = default;
 
    DbnBase& operator = (DbnBase&&) = default;
 
 
 
 
 
 
    void fill(const std::array<double, N> vals, const double weight = 1.0, const double fraction = 1.0) {
      fill(vals, weight, fraction, std::make_index_sequence<N>{});
    }
 
    template <typename... Args>
    void fill(Args&&... args) {
      // N + 2 - dimension + 2 default parameters
      static_assert(sizeof...(args) <= N + 2,
                    "Value's dimension doesn't match distribution's dimension.");
      std::array<double, sizeof...(args)> vals = {{args...}};
      double weight = 1.0;
      double fraction = 1.0;
      if (vals.size() > N) {
        weight = vals[N];
        if (vals.size() > N + 1)  fraction = vals[N + 1];
      }
      const double sf = fraction * weight;
      _numEntries += fraction;
      _sumW.at(0) += sf;
      _sumW2.at(0) += fraction*sqr(weight);
      for (unsigned int i = 0; i < N; ++i) {
        _sumW.at(i + 1) += sf * vals.at(i);
        _sumW2.at(i + 1) += sf * sqr(vals.at(i));
      }
 
      if constexpr(N >= 2) {
        size_t idx = 0;
        for (size_t i = 0; i < (N-1); ++i) {
          for (size_t j = i+1; j < N; ++j) {
            _sumWcross.at(idx++) += sf * vals.at(i) * vals.at(j);
          }
        }
      }
 
    }
 
 
    void set(const double numEntries, const std::array<double,N+1>& sumW,
                                      const std::array<double,N+1>& sumW2,
                                      const std::array<double,N*(N-1)/2>& sumWcross) {
      _numEntries = numEntries;
      _sumW = sumW;
      _sumW2 = sumW2;
      _sumWcross = sumWcross;
    }
 
    void set(const double numEntries, const std::vector<double>& sumW,
                                      const std::vector<double>& sumW2,
                                      const std::vector<double>& sumWcross = {}) {
      if (!(sumW.size() <= (N + 1) || sumW2.size() <= (N + 1) || sumWcross.size() <= (N*(N-1)/2)))
        throw UserError("Value's dimension doesn't match distribution's dimension.");
      _numEntries = numEntries;
      std::copy_n(std::make_move_iterator(sumW.begin()),      sumW.size(),      _sumW.begin());
      std::copy_n(std::make_move_iterator(sumW2.begin()),     sumW2.size(),     _sumW2.begin());
      std::copy_n(std::make_move_iterator(sumWcross.begin()), sumWcross.size(), _sumWcross.begin());
    }
 
 
    void set(const DbnBase<N>& other) {
      if (this != &other)  *this = other;
    }
 
    void set(DbnBase<N>&& other) {
      if (this != &other)  *this = std::move(other);
    }
 
    void reset() {
      _numEntries= 0;
      _sumW.fill(0);
      _sumW2.fill(0);
      _sumWcross.fill(0);
    }
 
 
    void scaleW(const double scalefactor) {
      _sumW.at(0) *= scalefactor;
      _sumW2.at(0) *= sqr(scalefactor);
      for (size_t i = 0; i< N; ++i) {
        // first- and second-order moments
        _sumW.at(i+1) *= scalefactor;
        _sumW2.at(i+1) *= scalefactor;
      }
      for (size_t i = 0; i < _sumWcross.size(); ++i) {
        _sumWcross.at(i) *= scalefactor;
      }
    }
 
    void scale(const size_t dim, const double factor) {
      if (dim >= N)
        throw RangeError("Dimension index should be less than "+std::to_string(N));
      _sumW.at(dim+1) *= factor;
      _sumW2.at(dim+1) *= sqr(factor);
      size_t idx = 0;
      for (size_t i = 0; i < (N-1); ++i) {
        for (size_t j = i+1; j < N; ++j) {
          if (i == dim || j == dim) {
            _sumWcross.at(idx++) *= factor;
          }
        }
      } // end of double loop
    }
 
 
 
    public:
 
 
 
 
    double errW() const;
 
    double relErrW() const;
 
    double mean(const size_t i) const;
 
    double variance(const size_t i) const;
 
    double stdDev(const size_t i) const { return std::sqrt(variance(i)); }
 
    double stdErr(const size_t i) const;
 
    double relStdErr(const size_t i) const;
 
    double RMS(const size_t i) const;
 
 
 
 
 
    double numEntries() const {
      return _numEntries;
    }
 
    double effNumEntries() const {
      if (_sumW2.at(0) == 0) return 0;
      return _sumW.at(0)*_sumW.at(0) / _sumW2.at(0);
    }
 
    double sumW() const {
      return _sumW.at(0);
    }
 
    double sumW2() const {
      return _sumW2.at(0);
    }
 
    double sumW(const size_t i) const {
      return _sumW.at(i);
    }
 
    double sumW2(const size_t i) const {
      return _sumW2.at(i);
    }
 
    template<size_t dim = N, typename = std::enable_if_t<(dim >= 2)>>
    double crossTerm(const size_t A1, const size_t A2) const {
      if (A1 >= N || A2 >= N)  throw RangeError("Invalid axis int, must be in range 0..dim-1");
      if (A1 >= A2)  throw RangeError("Indices need to be different for cross term");
 
      size_t idx = 0;
      for (size_t i = 0; i < (N-1); ++i) {
        for (size_t j = i+1; j < N; ++j) {
          if (i == A1 && j == A2)  break;
          ++idx;
        }
        if (i == A1)  break;
      }
 
      return _sumWcross.at(idx);
    }
 
    size_t dim() const { return N; }
 
 
 
    std::string toString() const {
      std::string res ="";
      res += ("numEntries="+ std::to_string(_numEntries))  ;
      if (sumW()) {
        res += (", Mean="+ std::to_string(mean(0)))  ;
        res += (", RMS="+ std::to_string(RMS(0)))  ;
      }
      return res;
    }
 
 
    // @{
 
    template<typename... size_t>
    DbnBase<sizeof...(size_t)> reduceTo(size_t... axes) const {
      if constexpr (N == sizeof...(axes)) { return *this; } // keep all axes
      else {
        // only select the axis moments to keep (0 is sumW, 1...N are the axis moments)
        std::array<double, sizeof...(axes)+1> newSumW = { _sumW[0], _sumW[axes+1]... };
        std::array<double, sizeof...(axes)+1> newSumW2 = { _sumW2[0], _sumW2[axes+1]... };
        if constexpr (sizeof...(axes) < 2) {
          return DbnBase<sizeof...(axes)>(_numEntries, newSumW, newSumW2);
        }
        else {
          constexpr auto newDim = sizeof...(axes);
          std::array<double, newDim*(newDim-1)/2> newSumWcross;
          unsigned int idx = 0;
          for (unsigned int i : {axes...}) {
            for (unsigned int j : {axes...}) {
              if (i < j) { newSumWcross.at(idx) = crossTerm(i, j); ++idx; }
            }
          }
          return DbnBase<sizeof...(axes)>(_numEntries, newSumW, newSumW2, newSumWcross);
        }
      }
    }
 
    template<size_t... Is>
    DbnBase<sizeof...(Is)> reduceTo(std::index_sequence<Is...>) const {
      return reduceTo(Is...);
    }
 
    template<size_t... axes>
    auto reduce() const {
      constexpr auto new_axes = AxisFilter(std::make_index_sequence<N>{},
        [](size_t i) { return !has_match(i, std::integer_sequence<size_t, axes...>{}); }
      );
      return reduceTo(new_axes);
    }
 
    // @}
 
 
 
    size_t _lengthContent() const noexcept {
      return DataSize::value;
    }
 
    std::vector<double> _serializeContent() const noexcept {
      std::vector<double> rtn;
      rtn.reserve(DataSize::value);
      rtn.insert(rtn.end(), _sumW.begin(),  _sumW.end());
      rtn.insert(rtn.end(), _sumW2.begin(), _sumW2.end());
      rtn.insert(rtn.end(), _sumWcross.begin(), _sumWcross.end());
      rtn.push_back(_numEntries);
      return rtn;
    }
 
    void _deserializeContent(const std::vector<double>& data) {
 
      if (data.size() != DataSize::value)
        throw UserError("Length of serialized data should be "+std::to_string(DataSize::value)+"!");
 
      auto itr = data.cbegin();
      std::copy_n(itr, _sumW.size(), _sumW.begin());
      std::copy_n(itr+N+1, _sumW2.size(), _sumW2.begin());
      std::copy_n(itr+2*(N+1), _sumWcross.size(), _sumWcross.begin());
      _numEntries = *(itr + 2*(N+1) + (N*(N-1)/2));
 
    }
 
    // @}
 
 
 
    DbnBase& operator += (const DbnBase& d);
 
    DbnBase& operator += (DbnBase&& d);
 
    DbnBase& operator -= (const DbnBase& d);
 
    DbnBase& operator -= (DbnBase&& d);
 
 
    private:
 
 
 
    double _numEntries;
 
    std::array<double, N+1> _sumW; //indexed by dimension
 
    std::array<double, N+1> _sumW2; //indexed by dimension
 
    std::array<double, N*(N-1)/2> _sumWcross;
 
 
 
 
    template <size_t... I>
    void fill(const std::array<double, N> vals, const double weight, const double fraction,
              std::index_sequence<I...>) {
      fill(vals[I]..., weight, fraction);
    }
 
  };
 
 
  template <size_t N>
  double DbnBase<N>::errW() const {
    return sqrt(sumW2(0));
  }
 
  template <size_t N>
  double DbnBase<N>::relErrW() const {
    if (effNumEntries() == 0) {
      return std::numeric_limits<double>::quiet_NaN();
    }
    return errW()/sumW(0);
  }
 
  template <size_t N>
  double DbnBase<N>::mean(const size_t i) const {
    return YODA::mean(sumW(i), sumW(0));
  }
 
 
  template <size_t N>
  double DbnBase<N>::variance(const size_t i) const {
    return YODA::variance(sumW(i), sumW(0), sumW2(i), sumW2(0));
  }
 
 
  template <size_t N>
  double DbnBase<N>::stdErr(const size_t i) const {
    return YODA::stdErr(sumW(i), sumW(0), sumW2(i), sumW2(0));
  }
 
 
  template <size_t N>
  double DbnBase<N>::relStdErr(const size_t i) const {
    if (effNumEntries() == 0) {
      return std::numeric_limits<double>::quiet_NaN();
    }
    return stdErr(i) / mean(i);
  }
 
 
  template <size_t N>
  double DbnBase<N>::RMS(const size_t i) const {
    return YODA::RMS(sumW2(i), sumW(0), sumW2());
  }
 
  template <size_t N>
  DbnBase<N>& DbnBase<N>::operator += (const DbnBase<N>& d) {
    _numEntries += d._numEntries;
    for (size_t i = 0; i <= N; ++i) {
      _sumW.at(i) += d._sumW.at(i);
      _sumW2.at(i) += d._sumW2.at(i);
    }
    for (size_t i = 0; i < _sumWcross.size(); ++i) {
      _sumWcross.at(i) += d._sumWcross.at(i);
    }
    return *this;
  }
 
 
  template <size_t N>
  DbnBase<N>& DbnBase<N>::operator += (DbnBase<N>&& d) {
    _numEntries += std::move(d._numEntries);
    for (size_t i = 0; i <= N; ++i) {
      _sumW.at(i) += std::move(d._sumW.at(i));
      _sumW2.at(i) += std::move(d._sumW2.at(i));
    }
    for (size_t i = 0; i < _sumWcross.size(); ++i) {
      _sumWcross.at(i) += std::move(d._sumWcross.at(i));
    }
    return *this;
  }
 
 
  template <size_t N>
  DbnBase<N>& DbnBase<N>::operator -= (const DbnBase<N>& d) {
    _numEntries -= d._numEntries;
    for (unsigned int i =0; i<= N ; ++i) {
      _sumW.at(i)  -= d._sumW.at(i);
      _sumW2.at(i) += d._sumW2.at(i);
    }
    for (size_t i = 0; i < _sumWcross.size(); ++i) {
      _sumWcross.at(i) -= d._sumWcross.at(i);
    }
    return *this;
  }
 
 
  template <size_t N>
  DbnBase<N>& DbnBase<N>::operator -= (DbnBase<N>&& d) {
    _numEntries -= std::move(d._numEntries);
    for (unsigned int i =0; i<= N ; ++i) {
      _sumW.at(i)  -= std::move(d._sumW.at(i));
      _sumW2.at(i) += std::move(d._sumW2.at(i));
    }
    for (size_t i = 0; i < _sumWcross.size(); ++i) {
      _sumWcross.at(i) -= std::move(d._sumWcross.at(i));
    }
    return *this;
  }
 
 
  template <size_t N>
  inline DbnBase<N> operator + (DbnBase<N> first, const DbnBase<N>& second) {
    first += second;
    return first;
  }
 
  template <size_t N>
  inline DbnBase<N> operator + (DbnBase<N> first, DbnBase<N>&& second) {
    first += std::move(second);
    return first;
  }
 
  template <size_t N>
  inline DbnBase<N> operator - (DbnBase<N> first, const DbnBase<N>& second) {
    first -= second;
    return first;
  }
 
  template <size_t N>
  inline DbnBase<N> operator - (DbnBase<N> first, DbnBase<N>&& second) {
    first -= std::move(second);
    return first;
  }
 
 
  template<size_t N>
  class Dbn : public DbnBase<N> {
    public:
    using BaseT = DbnBase<N>;
    using BaseT::BaseT;
    using BaseT::operator=;
 
  };
 
  template<>
  class Dbn<0> : public DbnBase<0> {
    public:
 
    using BaseT = DbnBase<0>;
    using BaseT::BaseT;
    using BaseT::operator=;
 
    Dbn() : BaseT() {}
 
    Dbn(const double numEntries, const double sumW, const double sumW2)
    : BaseT(numEntries, {sumW} , {sumW2}) { }
 
    Dbn(const BaseT& other) : BaseT(other) {}
 
    Dbn(BaseT&& other) : BaseT(std::move(other)) {}
 
    void fill(const std::array<double, 0>& vals, const double weight = 1.0, const double fraction = 1.0) {
      BaseT::fill(vals, weight, fraction);
    }
 
    void fill(const double weight=1.0, const double fraction=1.0) { BaseT::fill({}, weight, fraction); }
 
    void set(const double numEntries, const double sumW, const double sumW2) {
      BaseT::set(numEntries, {sumW}, {sumW2});
    }
 
  };
 
 
  template<>
  class Dbn<1> : public DbnBase<1>,
                 public XDbnMixin<Dbn<1>> {
    public:
 
    using BaseT = DbnBase<1>;
    using BaseT::BaseT;
    using BaseT::operator=;
 
    Dbn() : BaseT() {}
 
    Dbn(const double numEntries, const double sumW, const double sumW2, const double sumWX, const double sumWX2)
    : BaseT(numEntries, {sumW, sumWX} , {sumW2, sumWX2}) { }
 
    Dbn(const BaseT& other) : BaseT(other) {}
 
    Dbn(BaseT&& other) : BaseT(std::move(other)) {}
 
    void fill(const std::array<double, 1>& vals, const double weight = 1.0, const double fraction = 1.0) {
      BaseT::fill(vals, weight, fraction);
    }
 
    void fill(const double val, const double weight=1.0, const double fraction=1.0) {
      BaseT::fill({val}, weight, fraction);
    }
 
  };
 
 
  template<>
  class Dbn<2> : public DbnBase<2>,
                 public XDbnMixin<Dbn<2>>,
                 public YDbnMixin<Dbn<2>> {
    public:
 
    using BaseT = DbnBase<2>;
    using BaseT::BaseT;
    using BaseT::operator=;
 
    Dbn() : BaseT() {}
 
    Dbn(const double numEntries, const double sumW, const double sumW2, const double sumWX, const double sumWX2,
                                 const double sumWY, const double sumWY2, const double sumWXY)
    : BaseT(numEntries, {sumW, sumWX, sumWY} , {sumW2, sumWX2, sumWY2}, {sumWXY}) { }
 
    Dbn(const BaseT& other) : BaseT(other) {}
 
    Dbn(BaseT&& other) : BaseT(std::move(other)) {}
 
    void fill(const std::array<double, 2>& vals, const double weight = 1.0, const double fraction = 1.0) {
      BaseT::fill(vals, weight, fraction);
    }
 
    void fill(const double valX, const double valY, const double weight=1.0, const double fraction=1.0) {
      BaseT::fill({valX, valY}, weight, fraction);
    }
 
 
  };
 
 
  template<>
  class Dbn<3> : public DbnBase<3>,
                 public XDbnMixin<Dbn<3>>,
                 public YDbnMixin<Dbn<3>>,
                 public ZDbnMixin<Dbn<3>> {
    public:
 
    using BaseT = DbnBase<3>;
    using BaseT::BaseT;
    using BaseT::operator=;
 
    Dbn() : BaseT() {}
 
    Dbn(const double numEntries, const double sumW, const double sumW2,
                                 const double sumWX, const double sumWX2,
                                 const double sumWY, const double sumWY2,
                                 const double sumWZ, const double sumWZ2,
                                 const double sumWXY, const double sumWXZ, const double sumWYZ)
    : BaseT(numEntries, {sumW, sumWX, sumWY, sumWZ} , {sumW2, sumWX2, sumWY2, sumWZ2}, {sumWXY, sumWXZ, sumWYZ}) { }
 
    Dbn(const BaseT& other) : BaseT(other) {}
 
    Dbn(BaseT&& other) : BaseT(std::move(other)) {}
 
    void fill(const std::array<double, 3>& vals, const double weight = 1.0, const double fraction = 1.0) {
      BaseT::fill(vals, weight, fraction);
    }
 
    void fill(const double valX, const double valY, const double valZ, const double weight=1.0, const double fraction=1.0) {
      BaseT::fill({valX, valY, valZ}, weight, fraction);
    }
 
  };
 
  using Dbn0D = Dbn<0>;
  using Dbn1D = Dbn<1>;
  using Dbn2D = Dbn<2>;
  using Dbn3D = Dbn<3>;
}
 
#endif