Histogram.h

Go to the documentation of this file.

/***********************************************************************************\
* (c) Copyright 1998-2022 CERN for the benefit of the LHCb and ATLAS collaborations *
*                                                                                   *
* This software is distributed under the terms of the Apache version 2 licence,     *
* copied verbatim in the file "LICENSE".                                            *
*                                                                                   *
* In applying this licence, CERN does not waive the privileges and immunities       *
* granted to it by virtue of its status as an Intergovernmental Organization        *
* or submit itself to any jurisdiction.                                             *
\***********************************************************************************/
#pragma once
 
#include <Gaudi/Accumulators.h>
#include <Gaudi/MonitoringHub.h>
#include <GaudiKernel/HistoDef.h>
 
#include <array>
#include <cmath>
#include <fmt/format.h>
#include <nlohmann/json.hpp>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace Gaudi::Accumulators {
 
  namespace details {
    template <std::size_t, typename T>
    using alwaysT = T;
    // get a tuple of n types the given Type, or directly the type for n = 1
    template <typename Type, unsigned int ND>
    struct GetTuple;
    template <typename Type, unsigned int ND>
    using GetTuple_t = typename GetTuple<Type, ND>::type;
    template <typename Type, unsigned int ND>
    struct GetTuple {
      using type =
          decltype( std::tuple_cat( std::declval<std::tuple<Type>>(), std::declval<GetTuple_t<Type, ND - 1>>() ) );
    };
    template <typename Type>
    struct GetTuple<Type, 1> {
      using type = std::tuple<Type>;
    };
 
    inline void requireValidTitle( std::string_view sv ) {
      if ( !sv.empty() && ( std::isspace( sv.back() ) || std::isspace( sv.front() ) ) ) {
        throw GaudiException(
            fmt::format( "Histogram title \'{}\' has whitespace at front or back -- please remove", sv ),
            "Gaudi::Accumulators", StatusCode::FAILURE );
      }
    }
  } // namespace details
 
  struct ExtractWeight {
    template <typename Arithmetic>
    constexpr decltype( auto ) operator()( const std::pair<Arithmetic, Arithmetic>& v ) const noexcept {
      return v.second;
    }
  };
 
  struct WeightedProduct {
    template <typename Arithmetic>
    constexpr decltype( auto ) operator()( const std::pair<Arithmetic, Arithmetic>& v ) const noexcept {
      return v.first * v.second;
    }
  };
 
  struct WeightedSquare {
    template <typename Arithmetic>
    constexpr decltype( auto ) operator()( const std::pair<Arithmetic, Arithmetic>& v ) const noexcept {
      return v.first * v.first * v.second;
    }
  };
 
  template <atomicity Atomicity, typename Arithmetic>
  struct WeightedCountAccumulator
      : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, ExtractWeight> {
    using Base = GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, ExtractWeight>;
    using Base::Base;
    using Base::operator+=;
    WeightedCountAccumulator operator+=( const Arithmetic weight ) {
      *this += { Arithmetic{}, weight };
      return *this;
    }
    Arithmetic nEntries() const { return this->value(); }
  };
 
  template <atomicity Atomicity, typename Arithmetic>
  struct WeightedSumAccumulator
      : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, WeightedProduct> {
    using GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity,
                             WeightedProduct>::GenericAccumulator;
    Arithmetic sum() const { return this->value(); }
  };
 
  template <atomicity Atomicity, typename Arithmetic = double>
  struct WeightedSquareAccumulator
      : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, WeightedSquare> {
    using GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity,
                             WeightedSquare>::GenericAccumulator;
    Arithmetic sum2() const { return this->value(); };
  };
 
  template <atomicity Atomicity, typename Arithmetic>
  using WeightedAveragingAccumulator =
      AveragingAccumulatorBase<Atomicity, Arithmetic, WeightedCountAccumulator, WeightedSumAccumulator>;
 
  template <atomicity Atomicity, typename Arithmetic>
  using WeightedSigmaAccumulator =
      SigmaAccumulatorBase<Atomicity, Arithmetic, WeightedAveragingAccumulator, WeightedSquareAccumulator>;
 
  template <typename Arithmetic>
  struct Axis {
    Axis( unsigned int _nBins, Arithmetic _minValue, Arithmetic _maxValue, std::string _title = {},
          std::vector<std::string> _labels = {} )
        : nBins( _nBins )
        , minValue( _minValue )
        , maxValue( _maxValue )
        , title( std::move( _title ) )
        , labels( std::move( _labels ) )
        , ratio( _nBins / ( _maxValue - _minValue ) ) {
      details::requireValidTitle( title );
      for ( const auto& s : labels ) details::requireValidTitle( s );
    };
    explicit Axis( Gaudi::Histo1DDef const& def )
        : Axis( (unsigned int)def.bins(), def.lowEdge(), def.highEdge(), def.title() ){};
    unsigned int nBins;
    Arithmetic minValue, maxValue;
    std::string title;
    std::vector<std::string> labels;
    Arithmetic ratio;
 
    unsigned int index( Arithmetic value ) const {
      // In case we use integer as Arithmetic type, we cannot use ratio for computing indices,
      // as ratios < 1.0 will simply be 0, so we have to pay the division in such a case
      int idx;
      if constexpr ( std::is_integral_v<Arithmetic> ) {
        idx = ( ( value - minValue ) * nBins / ( maxValue - minValue ) ) + 1;
      } else {
        idx = std::floor( ( value - minValue ) * ratio ) + 1;
      }
      return idx < 0 ? 0 : ( (unsigned int)idx > nBins ? nBins + 1 : (unsigned int)idx );
    }
    bool inAcceptance( Arithmetic value ) const { return value >= minValue && value <= maxValue; }
  };
 
  template <typename Arithmetic>
  void to_json( nlohmann::json& j, const Axis<Arithmetic>& axis ) {
    j = nlohmann::json{ { "nBins", axis.nBins },
                        { "minValue", axis.minValue },
                        { "maxValue", axis.maxValue },
                        { "title", axis.title } };
    if ( !axis.labels.empty() ) { j["labels"] = axis.labels; }
  }
 
  template <typename Arithmetic, unsigned int ND, unsigned int NIndex = ND>
  struct HistoInputType : std::array<Arithmetic, ND> {
    // allow construction from set of values
    template <class... ARGS, typename = typename std::enable_if_t<( sizeof...( ARGS ) == NIndex )>>
    HistoInputType( ARGS... args ) : std::array<Arithmetic, ND>{ static_cast<Arithmetic>( args )... } {}
    // The change on NIndex == 1 allow to have simpler syntax in that case, that is no tuple of one item
    using ValueType          = HistoInputType<Arithmetic, NIndex == 1 ? 1 : ND, NIndex>;
    using AxisArithmeticType = Arithmetic;
    using InternalType       = std::array<Arithmetic, ND>;
    unsigned int computeIndex( const std::array<Axis<Arithmetic>, NIndex>& axis ) const {
      unsigned int index = 0;
      for ( unsigned int j = 0; j < NIndex; j++ ) {
        unsigned int dim = NIndex - j - 1;
        // compute local index for a given dimension
        int localIndex = axis[dim].index( ( *this )[dim] );
        // compute global index. Bins are stored in a column first manner
        index *= ( axis[dim].nBins + 2 );
        index += localIndex;
      }
      return index;
    }
    bool inAcceptance( const std::array<Axis<Arithmetic>, NIndex>& axis ) const {
      for ( unsigned int dim = 0; dim < NIndex; dim++ ) {
        if ( !axis[dim].inAcceptance( ( *this )[dim] ) ) return false;
      }
      return true;
    }
    auto forInternalCounter() { return 1ul; }
    template <typename AxisType, long unsigned NAxis>
    static unsigned int computeTotNBins( std::array<AxisType, NAxis> axis ) {
      unsigned int nTotBins = 1;
      for ( unsigned int i = 0; i < NAxis; i++ ) { nTotBins *= ( axis[i].nBins + 2 ); }
      return nTotBins;
    }
  };
 
  template <typename Arithmetic>
  class HistoInputType<Arithmetic, 1> {
  public:
    using ValueType          = HistoInputType;
    using AxisArithmeticType = Arithmetic;
    using InternalType       = Arithmetic;
    HistoInputType( Arithmetic a ) : value( a ) {}
    unsigned int computeIndex( const std::array<Axis<Arithmetic>, 1>& axis ) const { return axis[0].index( value ); }
    bool inAcceptance( const std::array<Axis<Arithmetic>, 1>& axis ) const { return axis[0].inAcceptance( value ); }
    Arithmetic& operator[]( int ) { return value; }
    operator Arithmetic() const { return value; }
    auto forInternalCounter() { return 1ul; }
    template <typename AxisType>
    static unsigned int computeTotNBins( std::array<AxisType, 1> axis ) {
      return axis[0].nBins + 2;
    }
 
  private:
    Arithmetic value;
  };
 
  template <typename Arithmetic, unsigned int ND, unsigned int NIndex = ND>
  struct WeightedHistoInputType : std::pair<HistoInputType<Arithmetic, ND, NIndex>, Arithmetic> {
    // The change on NIndex == 1 allow to have simpler syntax in that case, that is no tuple of one item
    using ValueType          = HistoInputType<Arithmetic, NIndex == 1 ? 1 : ND, NIndex>;
    using AxisArithmeticType = Arithmetic;
    using std::pair<HistoInputType<Arithmetic, ND, NIndex>, Arithmetic>::pair;
    unsigned int computeIndex( const std::array<Axis<Arithmetic>, NIndex>& axis ) const {
      return this->first.computeIndex( axis );
    }
    unsigned int inAcceptance( const std::array<Axis<Arithmetic>, NIndex>& axis ) const {
      return this->first.inAcceptance( axis );
    }
    auto forInternalCounter() { return std::pair( this->first.forInternalCounter(), this->second ); }
    template <typename AxisType, long unsigned NAxis>
    static unsigned int computeTotNBins( std::array<AxisType, NAxis> axis ) {
      return HistoInputType<Arithmetic, ND, NIndex>::computeTotNBins( axis );
    }
  };
 
  template <atomicity Atomicity, typename InputType, typename Arithmetic, typename ND,
            template <atomicity Ato, typename Arith> typename BaseAccumulatorT>
  class HistogramingAccumulatorInternal {
    template <atomicity, typename, typename, typename, template <atomicity, typename> typename>
    friend class HistogramingAccumulatorInternal;
 
  public:
    using BaseAccumulator    = BaseAccumulatorT<Atomicity, Arithmetic>;
    using AxisArithmeticType = typename InputType::AxisArithmeticType;
    using AxisType           = Axis<AxisArithmeticType>;
    template <std::size_t... Is>
    HistogramingAccumulatorInternal( details::GetTuple_t<AxisType, ND::value> axis, std::index_sequence<Is...> )
        : m_axis{ std::get<Is>( axis )... }
        , m_totNBins{ InputType::computeTotNBins( m_axis ) }
        , m_value( new BaseAccumulator[m_totNBins] ) {
      reset();
    }
    template <atomicity ato>
    HistogramingAccumulatorInternal(
        construct_empty_t,
        const HistogramingAccumulatorInternal<ato, InputType, Arithmetic, ND, BaseAccumulatorT>& other )
        : m_axis( other.m_axis )
        , m_totNBins{ InputType::computeTotNBins( m_axis ) }
        , m_value( new BaseAccumulator[m_totNBins] ) {
      reset();
    }
    [[deprecated( "Use `++h1[x]`, `++h2[{x,y}]`, etc. instead." )]] HistogramingAccumulatorInternal&
    operator+=( InputType v ) {
      accumulator( v.computeIndex( m_axis ) ) += v.forInternalCounter();
      return *this;
    }
    void reset() {
      for ( unsigned int index = 0; index < m_totNBins; index++ ) accumulator( index ).reset();
    }
    template <atomicity ato>
    void mergeAndReset( HistogramingAccumulatorInternal<ato, InputType, Arithmetic, ND, BaseAccumulatorT>& other ) {
      assert( m_totNBins == other.m_totNBins );
      for ( unsigned int index = 0; index < m_totNBins; index++ ) {
        accumulator( index ).mergeAndReset( other.accumulator( index ) );
      }
    }
    [[nodiscard]] auto operator[]( typename InputType::ValueType v ) {
      return Buffer<BaseAccumulatorT, Atomicity, Arithmetic>{ accumulator( v.computeIndex( m_axis ) ) };
    }
 
    auto& axis() const { return m_axis; }
    auto  nBins( unsigned int i ) const { return m_axis[i].nBins; }
    auto  minValue( unsigned int i ) const { return m_axis[i].minValue; }
    auto  maxValue( unsigned int i ) const { return m_axis[i].maxValue; }
    auto  binValue( unsigned int i ) const { return accumulator( i ).value(); }
    auto  nEntries( unsigned int i ) const { return accumulator( i ).nEntries(); }
    auto  totNBins() const { return m_totNBins; }
 
  private:
    BaseAccumulator& accumulator( unsigned int index ) const {
      assert( index < m_totNBins );
      return m_value[index];
    }
    std::array<AxisType, ND::value> m_axis;
    unsigned int m_totNBins;
    std::unique_ptr<BaseAccumulator[]> m_value;
  };
 
  template <atomicity Atomicity, typename Arithmetic, typename ND>
  using HistogramingAccumulator = HistogramingAccumulatorInternal<Atomicity, HistoInputType<Arithmetic, ND::value>,
                                                                  unsigned long, ND, IntegralAccumulator>;
 
  template <atomicity Atomicity, typename Arithmetic, typename ND>
  using WeightedHistogramingAccumulator =
      HistogramingAccumulatorInternal<Atomicity, WeightedHistoInputType<Arithmetic, ND::value>, Arithmetic, ND,
                                      WeightedCountAccumulator>;
 
  template <atomicity Atomicity, typename Arithmetic, typename ND>
  using ProfileHistogramingAccumulator =
      HistogramingAccumulatorInternal<Atomicity, HistoInputType<Arithmetic, ND::value + 1, ND::value>, Arithmetic, ND,
                                      SigmaAccumulator>;
 
  template <atomicity Atomicity, typename Arithmetic, typename ND>
  using WeightedProfileHistogramingAccumulator =
      HistogramingAccumulatorInternal<Atomicity, WeightedHistoInputType<Arithmetic, ND::value + 1, ND::value>,
                                      Arithmetic, ND, WeightedSigmaAccumulator>;
 
  template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type,
            template <atomicity, typename, typename> typename Accumulator, typename Seq>
  class HistogramingCounterBaseInternal;
  template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type,
            template <atomicity, typename, typename> typename Accumulator, std::size_t... NDs>
  class HistogramingCounterBaseInternal<ND, Atomicity, Arithmetic, Type, Accumulator, std::index_sequence<NDs...>>
      : public BufferableCounter<Atomicity, Accumulator, Arithmetic, std::integral_constant<int, ND>> {
  public:
    using Parent           = BufferableCounter<Atomicity, Accumulator, Arithmetic, std::integral_constant<int, ND>>;
    using AccumulatorType  = Accumulator<Atomicity, Arithmetic, std::integral_constant<int, ND>>;
    using NumberDimensions = std::integral_constant<unsigned int, ND>;
    inline static const std::string typeString{ std::string{ Type } + ':' + typeid( Arithmetic ).name() };
    template <typename OWNER>
    HistogramingCounterBaseInternal( OWNER* owner, std::string const& name, std::string const& title,
                                     details::GetTuple_t<typename AccumulatorType::AxisType, ND> axis )
        : Parent( owner, name, *this, axis, std::make_index_sequence<ND>{} ), m_title( title ) {
      details::requireValidTitle( m_title );
    }
    template <typename OWNER>
    HistogramingCounterBaseInternal( OWNER* owner, std::string const& name, std::string const& title,
                                     details::alwaysT<NDs, typename AccumulatorType::AxisType>... allAxis )
        : HistogramingCounterBaseInternal( owner, name, title, std::make_tuple( allAxis... ) ) {}
    using Parent::print;
    template <typename stream>
    stream& printImpl( stream& o, bool /*tableFormat*/ ) const {
      o << ND << "D Histogram with config ";
      for ( unsigned int i = 0; i < ND; i++ ) {
        o << this->nBins( i ) << " " << this->minValue( i ) << " " << this->maxValue( i ) << " ";
      }
      return o;
    }
    std::ostream& print( std::ostream& o, bool tableFormat = false ) const override {
      return printImpl( o, tableFormat );
    }
    MsgStream&  print( MsgStream& o, bool tableFormat = false ) const override { return printImpl( o, tableFormat ); }
    friend void reset( HistogramingCounterBaseInternal& c ) { c.reset(); }
    friend void mergeAndReset( HistogramingCounterBaseInternal& h, HistogramingCounterBaseInternal& o ) {
      h.mergeAndReset( o );
    }
    friend void to_json( nlohmann::json& j, HistogramingCounterBaseInternal const& h ) {
      // get all bin values and compute total nbEntries
      std::vector<typename AccumulatorType::BaseAccumulator::OutputType> bins;
      bins.reserve( h.totNBins() );
      unsigned long totNEntries{ 0 };
      for ( unsigned int i = 0; i < h.totNBins(); i++ ) {
        bins.push_back( h.binValue( i ) );
        totNEntries += h.nEntries( i );
      }
      // build json
      j = { { "type", std::string( Type ) + ":" + typeid( Arithmetic ).name() },
            { "title", h.m_title },
            { "dimension", ND },
            { "empty", totNEntries == 0 },
            { "nEntries", totNEntries },
            { "axis", h.axis() },
            { "bins", bins } };
    }
 
  private:
    std::string const m_title;
  };
  template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type,
            template <atomicity, typename, typename> typename Accumulator>
  using HistogramingCounterBase =
      HistogramingCounterBaseInternal<ND, Atomicity, Arithmetic, Type, Accumulator, std::make_index_sequence<ND>>;
 
  namespace naming {
    inline constexpr char histogramString[]                = "histogram:Histogram";
    inline constexpr char weightedHistogramString[]        = "histogram:WeightedHistogram";
    inline constexpr char profilehistogramString[]         = "histogram:ProfileHistogram";
    inline constexpr char weightedProfilehistogramString[] = "histogram:WeightedProfileHistogram";
  } // namespace naming
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using Histogram =
      HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString, HistogramingAccumulator>;
 
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using WeightedHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::weightedHistogramString,
                                                    WeightedHistogramingAccumulator>;
 
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using ProfileHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::profilehistogramString,
                                                   ProfileHistogramingAccumulator>;
 
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using WeightedProfileHistogram =
      HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::weightedProfilehistogramString,
                              WeightedProfileHistogramingAccumulator>;
 
} // namespace Gaudi::Accumulators