include/YODA/FillableStorage.h

0001 // -*- C++ -*-
0002 //
0003 // This file is part of YODA -- Yet more Objects for Data Analysis
0004 // Copyright (C) 2008-2024 The YODA collaboration (see AUTHORS for details)
0005 //
0006 #ifndef YODA_FillableStorage_H
0007 #define YODA_FillableStorage_H
0008
0009 #include "YODA/BinnedStorage.h"
0010 #include "YODA/Dbn.h"
0011
0012 namespace YODA {
0013
0014   /// @brief Type to adapter mapping used when user didn't provide type adapter.
0015   template <size_t FillDim, typename BinT>
0016   struct defaultAdapter; /// @note Error points here if there is no default adapter.
0017
0018
0019   /// @brief Default fill adapter for binned type double
0020   template <size_t FillDim, template<size_t, typename, typename> class BinT, typename BinningT, size_t N>
0021   struct defaultAdapter<FillDim, BinT<N, double, BinningT>> {
0022
0023     using AdapterT = std::function<void(
0024       BinT<N, double, BinningT>&,
0025       typename BinningT::EdgeTypesTuple&&, double, double)>;
0026
0027     AdapterT _adapter = [](auto& storedNumber, auto&& /* coords */,
0028                            double weight, double /* fraction */) {
0029       storedNumber = storedNumber + weight;
0030     };
0031   };
0032
0033   namespace {
0034     /// @brief Aux method to create a tuple of doubles,
0035     /// used as additional Dbn coordinates in profiles
0036     template <size_t... Is>
0037     constexpr auto dblPadding(std::index_sequence<Is...>) {
0038       return std::tuple< std::decay_t<decltype((void)Is, std::declval<double>())>...>();
0039     }
0040   }
0041
0042   /// @brief Default fill adapter for binned type Dbn<N>
0043   ///
0044   /// @note When non fundamental type is used, template parameter differ.
0045   template <size_t DbnN, template<size_t, typename, typename> class BinT, typename BinningT, size_t N>
0046   struct defaultAdapter<DbnN, BinT<N, Dbn<DbnN>, BinningT>> {
0047     static_assert(DbnN >= N, "Dimension of the Dbn needs to be at least as high as binning dimension!");
0048
0049     /// The fill coordinates are a tuple of
0050     /// -) the bin edge types (in case of histograms)
0051     /// -) the bin edge types and additional doubles (in case of profiles)
0052     using FillCoords = decltype(std::tuple_cat(std::declval<typename BinningT::EdgeTypesTuple>(),
0053                      dblPadding(std::make_index_sequence<DbnN-BinningT::Dimension::value>{})));
0054
0055     using AdapterT = std::function<void(BinT<N, Dbn<DbnN>, BinningT>&, FillCoords&&, double, double)>;
0056
0057     /// @note This adapter nullifies discrete coordinates and transforms
0058     /// coordinate tuple to std::array<double, N> to pass it to the
0059     /// stored Dbn object. The coordinate is moved to avoid copying.
0060     /// In cases when continuous axis' edges are not of double type
0061     /// that may positively affect performance.
0062     AdapterT _adapter = [](auto& dbn, auto&& coords, double weight, double fraction) {
0063       using CoordsArrT = std::array<double, DbnN>;
0064       CoordsArrT dblCoords;
0065       size_t binIdx = dbn.index();
0066       auto nullifyDiscrete = [&dblCoords, &binIdx, &coords](auto I){
0067         using isArithmetic = typename BinningT::template is_Arithmetic<I>;
0068         dblCoords[I] = nullifyIfDiscCoord(std::move(std::get<I>(coords)),
0069                                           std::integral_constant<bool,
0070                                           isArithmetic::value>{}, binIdx);
0071       };
0072       MetaUtils::staticFor<BinningT::Dimension::value>(nullifyDiscrete);
0073       if constexpr (DbnN > BinningT::Dimension::value)
0074         dblCoords[N] = std::move(std::get<N>(coords));
0075       dbn.fill(std::move(dblCoords), weight, fraction);
0076     };
0077   };
0078
0079
0080   /// @brief FillableStorage, introduces FillAdapterT on top of BinnedStorage base class
0081   ///
0082   /// @note The additional abstraction layer is necessary to distinguish between binned objects
0083   /// that are "live" (i.e. fillable/incrementable, like a Histo1D) and those that are "dead"
0084   /// (e.g. a binned set of cross-section measurement points).
0085   template <size_t FillDim, typename BinContentT, typename... AxisT>
0086   class FillableStorage
0087         : public BinnedStorage<BinContentT, AxisT...> {
0088
0089     static_assert((FillDim >= sizeof...(AxisT)),
0090                   "Fill dimension should be at least as large as the binning dimension.");
0091
0092   protected:
0093     /// @brief Convenience alias to be used in constructor
0094     using BaseT = BinnedStorage<BinContentT, AxisT...>;
0095     using BinningT = typename BaseT::BinningT;
0096     using BinT = Bin<sizeof...(AxisT), BinContentT, BinningT>;
0097     using AdapterWrapperT = defaultAdapter<FillDim, BinT>;
0098     using FillableT = FillableStorage<FillDim, BinContentT, AxisT...>;
0099     using FillCoordsT = decltype(std::tuple_cat(std::declval<typename BinningT::EdgeTypesTuple>(),
0100                                 dblPadding(std::make_index_sequence<FillDim-sizeof...(AxisT)>{})));
0101
0102   public:
0103
0104     /// @brief Type of the fill coordinates
0105     using FillType = FillCoordsT;
0106
0107     /// @brief Fill dimension
0108     using FillDimension = std::integral_constant<size_t, FillDim>;
0109
0110     /// @brief Adapter type (type of lambda used to access stored object).
0111     ///
0112     /// BinT is a stored object type;
0113     /// FillCoordsT is a coordinates init list (or tuple) type;
0114     /// Doubles are for weight and fraction correspondingly.
0115     using FillAdapterT = std::function<void(BinT&, FillCoordsT&&, double, double)>;
0116
0117     /// @name Constructors
0118     // @{
0119
0120     /// @brief Nullary constructor for unique pointers etc.
0121     FillableStorage(FillAdapterT adapter = AdapterWrapperT()._adapter)
0122         : BaseT(), _fillAdapter(adapter), _nancount(0), _nansumw(0.), _nansumw2(0.) { }
0123
0124     /// @brief Constructs FillableStorage from Binning.
0125     FillableStorage(const BinningT& binning, FillAdapterT adapter = AdapterWrapperT()._adapter)
0126         : BaseT(binning), _fillAdapter(adapter), _nancount(0), _nansumw(0.), _nansumw2(0.) { }
0127
0128     /// @brief Constructs FillableStorage from Binning. Rvalue.
0129     FillableStorage(BinningT&& binning, FillAdapterT adapter = AdapterWrapperT()._adapter)
0130         : BaseT(std::move(binning)), _fillAdapter(adapter), _nancount(0), _nansumw(0.), _nansumw2(0.) { }
0131
0132     /// @brief Constructs binning from an adapter and vectors of axes' edges
0133     FillableStorage(const std::vector<AxisT>&... edges, FillAdapterT adapter = AdapterWrapperT()._adapter)
0134         : BaseT(edges...), _fillAdapter(adapter), _nancount(0), _nansumw(0.), _nansumw2(0.) { }
0135
0136     /// @brief Constructs binning from an adapter and Rvalue vectors of axes' edges
0137     FillableStorage(std::vector<AxisT>&&... edges, FillAdapterT adapter = AdapterWrapperT()._adapter)
0138         : BaseT(std::move(edges)...), _fillAdapter(adapter), _nancount(0), _nansumw(0.), _nansumw2(0.) { }
0139
0140     /// @brief Constructs binning from an adapter and a sequence of axes
0141     FillableStorage(const Axis<AxisT>&... axes, FillAdapterT adapter = AdapterWrapperT()._adapter)
0142         : BaseT(axes...), _fillAdapter(adapter), _nancount(0), _nansumw(0.), _nansumw2(0.) { }
0143
0144     /// @brief Constructs binning from an adapter and a sequence of Rvalue axes
0145     FillableStorage(Axis<AxisT>&&... axes, FillAdapterT adapter = AdapterWrapperT()._adapter)
0146         : BaseT(std::move(axes)...), _fillAdapter(adapter), _nancount(0), _nansumw(0.), _nansumw2(0.) { }
0147
0148     /// @brief Copy constructor.
0149     FillableStorage(const FillableStorage& other)
0150         : BaseT(other), _fillAdapter(other._fillAdapter),
0151           _nancount(other._nancount), _nansumw(other._nansumw), _nansumw2(other._nansumw2) { }
0152
0153     /// @brief Move constructor.
0154     FillableStorage(FillableStorage&& other)
0155         : BaseT(std::move(other)), _fillAdapter(std::move(other._fillAdapter)),
0156           _nancount(std::move(other._nancount)), _nansumw(std::move(other._nansumw)),
0157           _nansumw2(std::move(other._nansumw2)) { }
0158
0159     // @}
0160
0161     /// @name Methods
0162     // @{
0163
0164     /// @brief Triggers fill adapter on the bin corresponding to coords
0165     ///
0166     /// @note Accepts coordinates only as rvalue tuple. The tuple members
0167     /// are then moved (bringing tuple member to unspecified state) later in adapters.
0168     template <size_t... Is>
0169     int fill(FillCoordsT&& coords, std::index_sequence<Is...>,
0170              const double weight = 1.0, const double fraction = 1.0) noexcept {
0171
0172       // make sure the user isn't trying to fill with NaN ...
0173       // include all fill coordinates here
0174       if (containsNan(coords)) {
0175         _nancount += 1;
0176         _nansumw += weight*fraction;
0177         _nansumw2 += sqr(weight*fraction);
0178         return -1;
0179       }
0180       // select binned coordinates (possibly a subset of fill coordinates)
0181       auto binCoords = std::tuple<AxisT...>(std::get<Is>(coords)...);
0182       const size_t binIdx = FillableT::_binning.globalIndexAt(binCoords);
0183       _fillAdapter(BaseT::bin(binIdx), std::move(coords), weight, fraction);
0184       return int(binIdx);
0185     }
0186
0187     /// @brief Triggers fill adapter on the bin corresponding to coords
0188     int fill(FillCoordsT&& coords, const double weight = 1.0, const double fraction = 1.0) noexcept {
0189       return fill(std::move(coords), std::make_index_sequence<sizeof...(AxisT)>{}, weight, fraction);
0190     }
0191
0192
0193     /// @name Utilities
0194     // @{
0195
0196     /// @brief Returns the dimension of the filling tuple
0197     size_t fillDim() const { return FillDim; }
0198
0199     size_t nanCount() const { return _nancount; }
0200
0201     double nanSumW() const { return _nansumw; }
0202
0203     double nanSumW2() const { return _nansumw2; }
0204
0205     void setNanLog(size_t count, double sumw, double sumw2) {
0206       _nancount = count;
0207       _nansumw  = sumw;
0208       _nansumw2 = sumw2;
0209     }
0210
0211     /// @brief Reset the Fillable.
0212     ///
0213     /// Keep the binning but set all bin contents and related quantities to zero
0214     void reset() noexcept {
0215       _nancount = 0;
0216       _nansumw = _nansumw2 = 0.;
0217       BaseT::clearBins();
0218     }
0219
0220     // @}
0221
0222     /// @name Operators
0223     // @{
0224
0225     /// @brief Copy assignment
0226     FillableStorage& operator = (const FillableStorage& other) noexcept {
0227       if (this != &other) {
0228         _fillAdapter = other._fillAdapter;
0229         _nancount    = other._nancount;
0230         _nansumw     = other._nansumw;
0231         _nansumw2    = other._nansumw2;
0232         BaseT::operator=(other);
0233       }
0234       return *this;
0235     }
0236
0237     /// @brief Move assignment
0238     FillableStorage& operator = (FillableStorage&& other) noexcept {
0239       if (this != &other) {
0240         _fillAdapter = std::move(other._fillAdapter);
0241         _nancount    = std::move(other._nancount);
0242         _nansumw     = std::move(other._nansumw);
0243         _nansumw2    = std::move(other._nansumw2);
0244         BaseT::operator=(std::move(other));
0245       }
0246       return *this;
0247     }
0248
0249     /// @brief Add another BinnedStorage to this one.
0250     FillableStorage& operator += (const FillableStorage& other) {
0251       if (*this != other)
0252         throw std::logic_error("YODA::BinnedStorage<" + std::to_string(sizeof...(AxisT)) +\
0253                                ">: Cannot add BinnedStorages with different binnings.");
0254       size_t i = 0;
0255       for (auto& bin : FillableT::bins(true)) {
0256         bin += other.bin(i++);
0257       }
0258
0259       return *this;
0260     }
0261
0262     /// @brief Subtract another BinnedStorage from this one.
0263     FillableStorage& operator -= (const FillableStorage& other) {
0264       if (*this != other)
0265         throw std::logic_error("YODA::FillableStorage<" + std::to_string(sizeof...(AxisT)) +\
0266                                ">: Cannot substract FillableStorages with different binnings.");
0267
0268       size_t i = 0;
0269       for (auto& bin : FillableT::bins(true)) {
0270         bin -= other.bin(i++);
0271       }
0272
0273       return *this;
0274     }
0275
0276     // @}
0277
0278     private:
0279
0280       /// @brief Adapter used to access stored objects
0281       FillAdapterT _fillAdapter;
0282
0283       size_t _nancount;
0284
0285       double _nansumw, _nansumw2;
0286
0287   };
0288
0289 } // namespace YODA
0290
0291 #endif