Accumulators.h

Go to the documentation of this file.

/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergei Zharko [committer] */
 
 
#pragma once
 
#include <boost/histogram/fwd.hpp>  // for weighted_mean<>
#include <boost/histogram/make_histogram.hpp>
#include <boost/serialization/access.hpp>
 
#include <cmath>
#include <type_traits>
 
namespace boost::histogram
{
  template<class ValueType>
  class RootStyleProfileAccumulator {
   public:
    using Value_t    = ValueType;
    using ConstRef_t = const Value_t&;
    using Weight_t   = weight_type<Value_t>;
 
    RootStyleProfileAccumulator() = default;
 
    template<class OtherValueType>
    RootStyleProfileAccumulator(const RootStyleProfileAccumulator<OtherValueType>& other)
      : fSumWV(other.fSumWV)
      , fSumWV2(other.fSumWV2)
      , fSumW(other.fSumW)
      , fSumW2(other.fSumW2)
    {
    }
 
    RootStyleProfileAccumulator(ConstRef_t sumWV, ConstRef_t sumWV2, ConstRef_t sumW, ConstRef_t sumW2)
      : fSumWV(sumWV)
      , fSumWV2(sumWV2)
      , fSumW(sumW)
      , fSumW2(sumW2)
    {
    }
 
    void operator()(ConstRef_t v) { operator()(weight(1.), v); }
 
    void operator()(const Weight_t& w, ConstRef_t v)
    {
      fSumWV += w.value * v;
      fSumWV2 += w.value * v * v;
      fSumW += w.value;
      fSumW2 += w.value * w.value;
    }
 
    RootStyleProfileAccumulator& operator+=(const RootStyleProfileAccumulator& rhs)
    {
      if (rhs.fSumW == 0) {
        return *this;
      }
      fSumWV += rhs.fSumWV;
      fSumWV2 += rhs.fSumWV2;
      fSumW += rhs.fSumW;
      fSumW2 += rhs.fSumW2;
      return *this;
    }
 
    RootStyleProfileAccumulator& operator*=(ConstRef_t s) noexcept
    {
      fSumWV *= s;
      fSumWV2 *= s * s;
      return *this;
    }
 
    bool operator==(const RootStyleProfileAccumulator& rhs) const noexcept
    {
      return fSumWV == rhs.fSumWV && fSumWV2 == rhs.fSumWV2 && fSumW == rhs.fSumW && fSumW2 == rhs.fSumW2;
    }
 
    bool operator!=(const RootStyleProfileAccumulator& rhs) const noexcept { return !operator==(rhs); }
 
    ConstRef_t GetSumW() const noexcept { return fSumW; }
 
    ConstRef_t GetSumW2() const noexcept { return fSumW2; }
 
    ConstRef_t GetSumWV() const noexcept { return fSumWV; }
 
    ConstRef_t GetSumWV2() const noexcept { return fSumWV2; }
 
    Value_t GetMean() const noexcept { return fSumW ? (fSumWV / fSumW) : 0; }
 
    Value_t GetEffCount() const noexcept { return fSumW2 ? (fSumW * fSumW / fSumW2) : 0; }
 
    Value_t GetVariance() const noexcept { return fSumW ? (fSumWV2 / fSumW - GetMean() * GetMean()) : 0; }
 
    Value_t GetStdDev() const noexcept { return std::sqrt(GetVariance()); }
 
    Value_t GetSEM() const noexcept { return fSumW ? GetStdDev() / std::sqrt(GetEffCount()) : 0; }
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& make_nvp("fSumWV", fSumWV);
      ar& make_nvp("fSumWV2", fSumWV2);
      ar& make_nvp("fSumW", fSumW);
      ar& make_nvp("fSumW2", fSumW2);
    }
 
    Value_t fSumWV{};   
    Value_t fSumWV2{};  
    Value_t fSumW{};    
    Value_t fSumW2{};   
  };
}  // namespace boost::histogram
 
#ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
namespace std
{
  template<class T, class U>
  struct common_type<boost::histogram::RootStyleProfileAccumulator<T>,
                     boost::histogram::RootStyleProfileAccumulator<U>> {
    using type = boost::histogram::RootStyleProfileAccumulator<common_type_t<T, U>>;
  };
}  // namespace std
#endif
 
namespace boost::histogram
{
  template<typename T>
  using RootStyleProfileStorage = dense_storage<RootStyleProfileAccumulator<T>>;
 
  template<class Axis, class... Axes, class = detail::requires_axis<Axis>>
  auto MakeRootStyleProfile(Axis&& axis, Axes&&... axes)
  {
    return make_histogram_with(RootStyleProfileStorage<double>(), std::forward<Axis>(axis),
                               std::forward<Axes>(axes)...);
  }
}  // namespace boost::histogram