Histogram.h

Go to the documentation of this file.

/***********************************************************************************\
* (c) Copyright 1998-2021 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 <ROOT/RHist.hxx>
 
#include <boost/format.hpp>
#include <boost/histogram/algorithm/sum.hpp>
#include <boost/histogram/fwd.hpp>
#include <gsl/span>
 
#include <nlohmann/json.hpp>
 
#include <array>
#include <cmath>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace {
  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>;
  };
} // namespace
 
namespace Gaudi::Accumulators {
 
  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 ) ){};
    unsigned int nBins;
    Arithmetic minValue, maxValue;
    std::string title;
    std::vector<std::string> labels;
    Arithmetic ratio;
 
    unsigned int index( Arithmetic value ) const {
      int idx = std::floor( ( value - minValue ) * ratio ) + 1;
      return idx < 0 ? 0 : ( (unsigned int)idx > nBins ? nBins + 1 : (unsigned int)idx );
    }
  };
 
  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>
    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;
    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;
    }
    auto forInternalCounter() { return this->operator[]( ND - 1 ); }
  };
 
  template <typename Arithmetic>
  class HistoInputType<Arithmetic, 1> {
  public:
    using ValueType          = HistoInputType;
    using AxisArithmeticType = Arithmetic;
    HistoInputType( Arithmetic a ) : value( a ) {}
    unsigned int computeIndex( const std::array<Axis<Arithmetic>, 1>& axis ) const { return axis[0].index( value ); }
    Arithmetic&  operator[]( int ) { return value; }
                 operator Arithmetic() const { return value; }
    auto         forInternalCounter() { return value; }
 
  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 );
    }
    auto forInternalCounter() { return std::pair( this->first.forInternalCounter(), this->second ); }
  };
 
  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;
    template <std::size_t... Is>
    HistogramingAccumulatorInternal( GetTuple_t<Axis<AxisArithmeticType>, ND::value> axis, std::index_sequence<Is...> )
        : m_axis{ std::get<Is>( axis )... }
        , m_totNBins{ computeTotNBins() }
        , 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{ computeTotNBins() }, m_value( new BaseAccumulator[m_totNBins] ) {
      reset();
    }
    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( std::move( other.accumulator( index ) ) );
      }
    }
    auto operator[]( typename InputType::ValueType v ) {
      return Buffer<BaseAccumulatorT, Atomicity, Arithmetic>{ accumulator( v.computeIndex( m_axis ) ) };
    }
 
  protected:
    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];
    }
    unsigned int computeTotNBins() const {
      unsigned int nTotBins = 1;
      for ( unsigned int i = 0; i < ND::value; i++ ) { nTotBins *= ( m_axis[i].nBins + 2 ); }
      return nTotBins;
    }
    std::array<Axis<AxisArithmeticType>, 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>>;
    template <typename OWNER>
    HistogramingCounterBaseInternal( OWNER* owner, std::string const& name, std::string const& title,
                                     GetTuple_t<Axis<Arithmetic>, ND> axis )
        : Parent( owner, name, std::string( Type ) + ":" + typeid( Arithmetic ).name(), axis,
                  std::make_index_sequence<ND>{} )
        , m_title( title ) {}
    template <typename OWNER>
    HistogramingCounterBaseInternal( OWNER* owner, std::string const& name, std::string const& title,
                                     alwaysT<NDs, Axis<Arithmetic>>... 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 ); }
    nlohmann::json toJSON() const override {
      // get all bin values and compute total nbEntries
      using Acc = Accumulator<Atomicity, Arithmetic, std::integral_constant<int, ND>>;
      std::vector<typename Acc::BaseAccumulator::OutputType> bins;
      bins.reserve( this->totNBins() );
      unsigned long totNEntries{ 0 };
      for ( unsigned int i = 0; i < this->totNBins(); i++ ) {
        bins.push_back( this->binValue( i ) );
        totNEntries += this->nEntries( i );
      }
      // build json
      return { { "type", std::string( Type ) + ":" + typeid( Arithmetic ).name() },
               { "title", m_title },
               { "dimension", ND },
               { "empty", totNEntries == 0 },
               { "nEntries", totNEntries },
               { "axis", this->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 {
    static const char histogramString[]                = "histogram:Histogram";
    static const char weightedHistogramString[]        = "histogram:WeightedHistogram";
    static const char profilehistogramString[]         = "histogram:ProfileHistogram";
    static const char weightedProfilehistogramString[] = "histogram:WeightedProfileHistogram";
  } // namespace
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using Histogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, histogramString, HistogramingAccumulator>;
 
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using WeightedHistogram =
      HistogramingCounterBase<ND, Atomicity, Arithmetic, weightedHistogramString, WeightedHistogramingAccumulator>;
 
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using ProfileHistogram =
      HistogramingCounterBase<ND, Atomicity, Arithmetic, profilehistogramString, ProfileHistogramingAccumulator>;
 
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using WeightedProfileHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, weightedProfilehistogramString,
                                                           WeightedProfileHistogramingAccumulator>;
 
} // namespace Gaudi::Accumulators