Histogram.h

Go to the documentation of this file.

/***********************************************************************************\
* (c) Copyright 1998-2019 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 <type_traits>
#include <utility>
 
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, const std::string& _title = "" )
        : nBins( _nBins )
        , minValue( _minValue )
        , maxValue( _maxValue )
        , title( _title )
        , ratio( _nBins / ( _maxValue - _minValue ) ){};
    unsigned int nBins;
    Arithmetic minValue, maxValue;
    std::string title;
    Arithmetic ratio;
  };
 
  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}};
  }
 
  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 dim = 0; dim < NIndex; dim++ ) {
        // compute local index for a given dimension
        int localIndex = std::floor( ( ( *this )[dim] - axis[dim].minValue ) * axis[dim].ratio ) + 1;
        localIndex     = ( localIndex < 0 )
                         ? 0
                         : ( ( (unsigned int)localIndex > axis[dim].nBins ) ? axis[dim].nBins + 1
                                                                            : (unsigned int)( localIndex ) );
        // compute global index. Bins are stored in a row first manner
        index = ( dim > 0 ? ( axis[dim - 1].nBins + 2 ) : 0 ) * 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 {
      int index = std::floor( ( value - axis[0].minValue ) * axis[0].ratio ) + 1;
      return index < 0 ? 0 : ( (unsigned int)index > axis[0].nBins ? axis[0].nBins + 1 : (unsigned int)index );
    }
    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( std::initializer_list<Axis<AxisArithmeticType>> axis, std::index_sequence<Is...> )
        : m_axis{{*( axis.begin() + Is )...}}
        , 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  ratio( unsigned int i ) const { return m_axis[i].ratio; }
    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>
  class HistogramingCounterBase
      : 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>
    HistogramingCounterBase( OWNER* owner, std::string const& name, std::string const& title,
                             std::initializer_list<Axis<Arithmetic>> axis )
        : Parent( owner, name, std::string( Type ) + ":" + typeid( Arithmetic ).name(), axis,
                  std::make_index_sequence<ND>{} )
        , m_title( title ) {}
    template <typename OWNER>
    HistogramingCounterBase( OWNER* owner, std::string const& name, std::string const& title, Axis<Arithmetic> axis )
        : HistogramingCounterBase( owner, name, title, {axis} ) {}
    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 int 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 == Arithmetic{}},
              {"nEntries", totNEntries},
              {"axis", this->axis()},
              {"bins", bins}};
    }
 
  private:
    std::string const m_title;
  };
 
  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