include/YODA/BinnedEstimate.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_BinnedEstimate_h
0007 #define YODA_BinnedEstimate_h
0008
0009 #include "YODA/AnalysisObject.h"
0010 #include "YODA/BinnedStorage.h"
0011 #include "YODA/Estimate.h"
0012 #include "YODA/Scatter.h"
0013 #include "YODA/Transformation.h"
0014 #include "YODA/Utils/BinnedUtils.h"
0015 #include <regex>
0016
0017 namespace YODA {
0018
0019   /// @note Anonymous namespace to limit visibility to this file
0020   namespace {
0021     /// @brief Nullify bin uncertainty for discrete axes
0022     template<size_t axisN, typename BinT>
0023     std::pair<double, double> nullifyIfDisc(const BinT& /* b */, const double /* val */, std::false_type, const double null = 0.0) {
0024       return { null, null };
0025     }
0026
0027     /// @brief Work out bin uncertainty with respect to
0028     /// central value for continuous axes
0029     template<size_t axisN, typename BinT>
0030     std::pair<double, double> nullifyIfDisc(const BinT& b, const double val, std::true_type, const double /* null */ = 0.0) {
0031       return { val - b.template min<axisN>(), b.template max<axisN>() - val };
0032     }
0033
0034     /// @brief Unify bin interval for discrete axes
0035     template<size_t axisN, typename BinT>
0036     double unifyIfDisc(const BinT& /* b */, std::false_type, const double unity = 1.0) {
0037       return unity;
0038     }
0039
0040     /// @brief Get bin width for continuous axes
0041     template<size_t axisN, typename BinT>
0042     double unifyIfDisc(const BinT& b, std::true_type, const double /* unity */ = 1.0) {
0043       return b.template width<axisN>();
0044     }
0045
0046     /// @brief Return bin index for non-integral discrete axes
0047     template<size_t axisN, typename BinT>
0048     double coordPicker(const BinT& b, std::false_type, std::false_type) {
0049       return b.index();
0050     }
0051
0052     /// @brief Return bin edge for integral discrete axes
0053     template<size_t axisN, typename BinT>
0054     double coordPicker(const BinT& b, std::true_type, std::false_type) {
0055       return b.template edge<axisN>();
0056     }
0057     /// @brief Return bin mid point for continuous axes
0058     template<size_t axisN, typename BinT>
0059     double coordPicker(const BinT& b, std::true_type, std::true_type) {
0060       return b.template mid<axisN>();
0061     }
0062   }
0063
0064
0065   /// Forward declaration
0066   template <typename... AxisT>
0067   class BinnedEstimate;
0068
0069   /// @brief EstimateStorage convenience class based on BinnedStorage.
0070   ///
0071   /// @note We use this abstraction layer to implement the bulk once and only once.
0072   /// The user-facing BinnedEstimate type will inherit all its methods from this
0073   /// base class along with a few axis-specifc mixins.
0074   template<typename... AxisT>
0075   class EstimateStorage
0076       : public BinnedStorage<Estimate, AxisT...>,
0077         public AnalysisObject {
0078   protected:
0079
0080     using BaseT = BinnedStorage<Estimate, AxisT...>;
0081     using BinningT = typename BaseT::BinningT;
0082     using BinT = typename BaseT::BinT;
0083
0084   public:
0085
0086     using BinType = BinT;
0087     using AnalysisObject::operator =;
0088
0089     /// @name Constructors
0090     /// @{
0091
0092     /// @brief Nullary constructor for unique pointers etc.
0093     EstimateStorage(const std::string& path = "", const std::string& title = "")
0094         : BaseT(), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0095
0096     /// @brief Constructor giving explicit bin edges as rvalue reference.
0097     ///
0098     /// Discrete axes have as many edges as bins.
0099     /// Continuous axes have bins.size()+1 edges, the last one
0100     /// being the upper bound of the last bin.
0101     EstimateStorage(std::vector<AxisT>&&... binsEdges, const std::string& path = "", const std::string& title = "")
0102         : BaseT(Axis<AxisT>(std::move(binsEdges))...), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0103
0104     /// @brief Constructor giving explicit bin edges as lvalue reference.
0105     ///
0106     /// Discrete axes have as many edges as bins.
0107     /// Continuous axes have bins.size()+1 edges, the last one
0108     /// being the upper bound of the last bin.
0109     EstimateStorage(const std::vector<AxisT>&... binsEdges, const std::string& path = "", const std::string& title = "")
0110         : BaseT(Axis<AxisT>(binsEdges)...), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0111
0112     /// @brief Constructor giving explicit bin edges as initializer list.
0113     ///
0114     /// Discrete axes have as many edges as bins.
0115     /// Continuous axes have bins.size()+1 edges, the last one
0116     /// being the upper bound of the last bin.
0117     EstimateStorage(std::initializer_list<AxisT>&&... binsEdges, const std::string& path = "", const std::string& title = "")
0118         : BaseT(Axis<AxisT>(std::move(binsEdges))...), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0119
0120     /// @brief Constructor giving range and number of bins.
0121     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT...>>
0122     EstimateStorage(const std::vector<size_t>& nBins,
0123                     const std::vector<std::pair<EdgeT, EdgeT>>& limitsLowUp,
0124                     const std::string& path = "", const std::string& title = "")
0125         //: BaseT( BinningT(nBins, limitsLowUp, std::make_index_sequence<sizeof...(AxisT)>{}) ),
0126         : BaseT( _mkBinning(nBins, limitsLowUp, std::make_index_sequence<sizeof...(AxisT)>{}) ),
0127           AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0128
0129     /// @brief Constructor given a binning type
0130     EstimateStorage(const BinningT& binning, const std::string& path = "", const std::string& title = "")
0131         : BaseT(binning), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0132
0133     /// @brief Constructor given an rvalue BinningT
0134     EstimateStorage(BinningT&& binning, const std::string& path = "", const std::string& title = "")
0135         : BaseT(std::move(binning)), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0136
0137     /// @brief Constructor given a scatter
0138     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT...>>
0139     EstimateStorage(const ScatterND<sizeof...(AxisT)+1>& s, const std::string& path = "", const std::string& title = "")
0140          : BaseT(_mkBinning(s, std::make_index_sequence<sizeof...(AxisT)>{})),
0141            AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0142
0143     /// @brief Copy constructor
0144     ///
0145     /// @todo Also allow title setting from the constructor?
0146     EstimateStorage(const EstimateStorage& other, const std::string& path = "") : BaseT(other),
0147         AnalysisObject(mkTypeString<AxisT...>(), path!=""? path : other.path(), other, other.title()) { }
0148
0149     /// @brief Move constructor
0150     ///
0151     /// @todo Also allow title setting from the constructor?
0152     EstimateStorage(EstimateStorage&& other, const std::string& path = "") : BaseT(std::move(other)),
0153         AnalysisObject(mkTypeString<AxisT...>(), path!=""? path : other.path(), other, other.title()) { }
0154
0155     /// @brief Make a copy on the stack
0156     EstimateStorage clone() const noexcept {
0157       return EstimateStorage(*this);
0158     }
0159
0160     /// @brief Make a copy on the heap
0161     EstimateStorage* newclone() const noexcept {
0162       return new EstimateStorage(*this);
0163     }
0164
0165     /// @}
0166
0167
0168     /// @name I/O
0169     /// @{
0170
0171   private:
0172
0173     // @brief Render information about this AO (private implementation)
0174     template<size_t... Is>
0175     void _renderYODA_aux(std::ostream& os, const int width, std::index_sequence<Is...>) const noexcept {
0176
0177       // print bin edges
0178       BaseT::_binning._renderYODA(os);
0179
0180       // Assemble union of error sources, as it's not guaranteed
0181       // that every bin has the same error breakdown
0182       const std::vector<std::string> labels = this->sources();
0183       if (labels.size()) {
0184         os << "ErrorLabels: [";
0185         for (size_t i = 0; i < labels.size(); ++i) {
0186           const std::string& src = labels[i];
0187           if (i)  os << ", ";
0188           os << std::quoted(src);
0189         }
0190         os << "]\n";
0191       }
0192
0193       // column header: content types
0194       os << std::setw(width) << std::left << "# value" << "\t";
0195       const int errwidth = std::max(int(std::to_string(labels.size()).size()+7), width); // "errDn(" + src + ")"
0196       for (size_t i = 0; i < labels.size(); ++i) {
0197         const std::string& src = labels[i];
0198         if (src.empty()) {
0199           os << std::setw(errwidth) << std::left << "totalDn" << "\t"
0200              << std::setw(errwidth) << std::left << "totalUp" << "\t";
0201         }
0202         else {
0203           os << std::setw(errwidth) << std::left << ("errDn(" + std::to_string(i+1) + ")") << "\t"
0204              << std::setw(errwidth) << std::left << ("errUp(" + std::to_string(i+1) + ")") << "\t";
0205         }
0206       }
0207       os << "\n";
0208
0209       for (const auto& b : BaseT::bins(true, true)) {
0210         os << std::setw(width) << std::left << b.val() << "\t"; // print value
0211         // print systs if available
0212         for (const std::string& src : labels) {
0213           if (!b.hasSource(src)) {
0214             os << std::setw(errwidth) << std::left << "---" << "\t"
0215                << std::setw(errwidth) << std::left << "---" << "\t";
0216             continue;
0217           }
0218           const auto& err = b.err(src);
0219           os << std::setw(errwidth) << std::left << err.first << "\t"
0220              << std::setw(errwidth) << std::left << err.second << "\t";
0221         }
0222         os << "\n";
0223       }
0224     }
0225
0226   public:
0227
0228     /// @name I/O
0229     /// @{
0230
0231     // @brief Render information about this AO
0232     void _renderYODA(std::ostream& os, const int width = 13) const noexcept {
0233       _renderYODA_aux(os, width, std::make_index_sequence<sizeof...(AxisT)>{});
0234     }
0235
0236     // @brief Render scatter-like information about this AO
0237     void _renderFLAT(std::ostream& os, const int width = 13) const noexcept {
0238       const ScatterND<sizeof...(AxisT)+1> tmp = mkScatter();
0239       tmp._renderYODA(os, width);
0240     }
0241
0242     /// @}
0243
0244
0245     /// @name Transformations
0246     /// @{
0247
0248     /// @brief Rescale as if all fill weights had been different by factor @a scalefactor.
0249     void scale(const double scalefactor) noexcept {
0250       setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
0251       for (auto& bin : BaseT::_bins) {
0252         bin.scale(scalefactor);
0253       }
0254     }
0255
0256
0257     /// @brief Merge every group of @a n bins, from start to end inclusive
0258     ///
0259     /// If the number of bins is not a multiple of @a n, the last @a m < @a n
0260     /// bins on the RHS will also be merged, as the closest possible approach to
0261     /// factor @n rebinning everywhere.
0262     ///
0263     /// @note Only visible bins are being rebinned. Underflow (index = 0) and
0264     /// overflow (index = numBins() + 1) are not included.
0265     template <size_t axisN>
0266     void rebinBy(unsigned int n, size_t begin=1, size_t end=UINT_MAX) {
0267       if (n < 1)  throw UserError("Rebinning requested in groups of 0!");
0268       if (!begin) throw UserError("Visible bins start with index 1!");
0269       if (end > BaseT::numBinsAt(axisN)+1)  end = BaseT::numBinsAt(axisN) + 1;
0270       for (size_t m = begin; m < end; ++m) {
0271         if (m > BaseT::numBinsAt(axisN)+1) break; // nothing to be done
0272         const size_t myend = (m+n-1 < BaseT::numBinsAt(axisN)+1) ? m+n-1 : BaseT::numBinsAt(axisN);
0273         if (myend > m) {
0274           BaseT::template mergeBins<axisN>({m, myend});
0275           end -= myend-m; //< reduce upper index by the number of removed bins
0276         }
0277       }
0278     }
0279
0280     /// @brief Overloaded alias for rebinBy
0281     template <size_t axisN>
0282     void rebin(unsigned int n, size_t begin=1, size_t end=UINT_MAX) {
0283       rebinBy<axisN>(n, begin, end);
0284     }
0285
0286     /// @brief Rebin to the given list of bin edges
0287     template <size_t axisN>
0288     void rebinTo(const std::vector<typename BinningT::template getAxisT<axisN>::EdgeT>& newedges) {
0289       if (newedges.size() < 2)
0290         throw UserError("Requested rebinning to an edge list which defines no bins");
0291       using thisAxisT = typename BinningT::template getAxisT<axisN>;
0292       using thisEdgeT = typename thisAxisT::EdgeT;
0293       // get list of shared edges
0294       thisAxisT& oldAxis = BaseT::_binning.template axis<axisN>();
0295       const thisAxisT newAxis(std::move(newedges));
0296       const std::vector<thisEdgeT> eshared = oldAxis.sharedEdges(newAxis);
0297       if (eshared.size() != newAxis.edges().size())
0298         throw BinningError("Requested rebinning to incompatible edges");
0299       // loop over new lower bin edges (= first bin index of merge range)
0300       for (size_t begin = 0; begin < eshared.size() - 1; ++begin) {
0301         // find index of upper edge along old axis
0302         // (subtracting 1 gives index of last bin to be merged)
0303         size_t end = oldAxis.index(eshared[begin+1]) - 1;
0304         // if the current edge is the last visible edge before the overflow
0305         // merge the remaining bins into the overflow
0306         if (begin == newAxis.numBins(true)-1)  end = oldAxis.numBins(true)-1;
0307         // merge this range
0308         if (end > begin)  BaseT::template mergeBins<axisN>({begin, end});
0309         if (eshared.size() == oldAxis.edges().size())  break; // we're done
0310       }
0311     }
0312
0313     /// @brief Overloaded alias for rebinTo
0314     template <size_t axisN>
0315     void rebin(const std::vector<typename BinningT::template getAxisT<axisN>::EdgeT>& newedges) {
0316       rebinTo<axisN>(std::move(newedges));
0317     }
0318
0319     /// @brief Reset the EstimateStorage
0320     ///
0321     /// Keep the binning but set all bin contents and related quantities to zero
0322     void reset() noexcept { BaseT::clearBins(); }
0323
0324     /// @brief Copy assignment
0325     EstimateStorage& operator = (const EstimateStorage& est) noexcept {
0326       if (this != &est) {
0327         AnalysisObject::operator = (est);
0328         BaseT::operator = (est);
0329       }
0330       return *this;
0331     }
0332
0333     /// Move assignment
0334     EstimateStorage& operator = (EstimateStorage&& est) noexcept {
0335       if (this != &est) {
0336         AnalysisObject::operator = (est);
0337         BaseT::operator = (std::move(est));
0338       }
0339       return *this;
0340     }
0341
0342     /// @brief Add two EstimateStorages
0343     ///
0344     /// @note Adding EstimateStorages will unset any ScaledBy
0345     /// attribute from previous calls to scale or normalize.
0346     EstimateStorage& add(const EstimateStorage& est,
0347                          const std::string& pat_uncorr="^stat|^uncor" ) {
0348       if (*this != est)
0349         throw BinningError("Arithmetic operation requires compatible binning!");
0350       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0351       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0352         BaseT::bin(i).add(est.bin(i), pat_uncorr);
0353       }
0354       BaseT::maskBins(est.maskedBins(), true);
0355       return *this;
0356     }
0357     //
0358     EstimateStorage& operator += (const EstimateStorage& est) {
0359       return add(est);
0360     }
0361
0362     /// @brief Add two (rvalue) EstimateStorages
0363     ///
0364     /// @note Adding EstimateStorages will unset any ScaledBy
0365     /// attribute from previous calls to scale or normalize.
0366     EstimateStorage& add(EstimateStorage&& est,
0367                          const std::string& pat_uncorr="^stat|^uncor" ) {
0368       if (*this != est)
0369         throw BinningError("Arithmetic operation requires compatible binning!");
0370       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0371       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0372         BaseT::bin(i).add(std::move(est.bin(i)), pat_uncorr);
0373       }
0374       BaseT::maskBins(est.maskedBins(), true);
0375       return *this;
0376     }
0377     //
0378     EstimateStorage& operator += (EstimateStorage&& est) {
0379       return add(std::move(est));
0380     }
0381
0382
0383     /// @brief Subtract one EstimateStorages from another one
0384     ///
0385     /// @note Subtracting EstimateStorages will unset any ScaledBy
0386     /// attribute from previous calls to scale or normalize.
0387     EstimateStorage& subtract(const EstimateStorage& est,
0388                               const std::string& pat_uncorr="^stat|^uncor" ) {
0389       if (*this != est)
0390         throw BinningError("Arithmetic operation requires compatible binning!");
0391       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0392       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0393         BaseT::bin(i).subtract(est.bin(i), pat_uncorr);
0394       }
0395       BaseT::maskBins(est.maskedBins(), true);
0396       return *this;
0397     }
0398     //
0399     EstimateStorage& operator -= (const EstimateStorage& est) {
0400       return subtract(est);
0401     }
0402
0403     /// @brief Subtract one (rvalue) EstimateStorages from another one
0404     EstimateStorage& subtract(EstimateStorage&& est,
0405                               const std::string& pat_uncorr="^stat|^uncor" ) {
0406       if (*this != est)
0407         throw BinningError("Arithmetic operation requires compatible binning!");
0408       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0409       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0410         BaseT::bin(i) -= std::move(est.bin(i), pat_uncorr);
0411       }
0412       BaseT::maskBins(est.maskedBins(), true);
0413       return *this;
0414     }
0415     //
0416     EstimateStorage& operator -= (EstimateStorage&& est) {
0417       return subtract(std::move(est));
0418     }
0419
0420     /// @}
0421
0422
0423     /// @name Binning info.
0424     /// @{
0425
0426     /// @brief Total dimension of the object (number of axes + estimate)
0427     size_t dim() const noexcept { return sizeof...(AxisT) + 1; }
0428
0429     /// @brief Returns the axis configuration
0430     std::string _config() const noexcept { return mkAxisConfig<AxisT...>(); }
0431
0432     /// @brief Templated version to get edges of axis N by value.
0433     /// +-inf edges are included.
0434     ///
0435     /// @note Needed by axis-specific verison from AxisMixin
0436     template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
0437     std::vector<E> edges(const bool includeOverflows = false) const noexcept {
0438       return BaseT::_binning.template edges<I>(includeOverflows);
0439     }
0440
0441     /// @brief Templated version to get bin widths of axis N by value.
0442     ///
0443     /// Overflows are included depending on @a includeOverflows
0444     /// Needed by axis-specific version from AxisMixin
0445     ///
0446     /// @note This version is only supported for continuous axes.
0447     template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
0448     std::enable_if_t<std::is_floating_point<E>::value, std::vector<E>>
0449     widths(const bool includeOverflows = false) const noexcept {
0450       return BaseT::_binning.template widths<I>(includeOverflows);
0451     }
0452
0453     /// @brief Get the lowest non-overflow edge of the axis
0454     ///
0455     /// @note This method is only supported for continuous axes
0456     template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
0457     std::enable_if_t<std::is_floating_point<E>::value, E> min() const noexcept {
0458       return BaseT::_binning.template min<I>();
0459     }
0460
0461     /// @brief Get the highest non-overflow edge of the axis
0462     ///
0463     /// @note This method is only supported for continuous axes
0464     template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
0465     std::enable_if_t<std::is_floating_point<E>::value, E> max() const noexcept {
0466       return BaseT::_binning.template max<I>();
0467     }
0468
0469     /// @}
0470
0471
0472     /// @name Whole EstimateStorage data
0473     /// @{
0474
0475     /// @brief Get list of central values
0476     std::vector<double> vals(const bool includeOverflows=false,
0477                              const bool includeMaskedBins=false) const {
0478       std::vector<double> rtn;
0479       rtn.reserve(BaseT::numBins(includeOverflows, includeMaskedBins));
0480       for (const auto& b : BaseT::bins(includeOverflows, includeMaskedBins)) {
0481         rtn.push_back(b.val());
0482       }
0483       return rtn;
0484     }
0485
0486     /// @brief Get list of error sources
0487     std::vector<std::string> sources() const {
0488       // Assemble union of error sources, as it's not guaranteed
0489       // that every bin has the same error breakdown
0490       std::vector<std::string> rtn;
0491
0492       for (const auto& b : BaseT::bins(true)) {
0493         std::vector<std::string> keys = b.sources();
0494         rtn.insert(std::end(rtn),
0495                    std::make_move_iterator(std::begin(keys)),
0496                    std::make_move_iterator(std::end(keys)));
0497       }
0498       std::sort(rtn.begin(), rtn.end());
0499       rtn.erase( std::unique(rtn.begin(), rtn.end()), rtn.end() );
0500
0501       return rtn;
0502     }
0503
0504     /// @brief Calculate the volume underneath the EstimateStorage
0505     ///
0506     /// @note overflow bins have infinite bin width
0507     double areaUnderCurve(const bool includeBinVol=true,
0508                           const bool includeOverflows=false,
0509                           const bool includeMaskedBins=false) const {
0510       double ret = 0.;
0511       for (const auto& b : BaseT::bins(includeOverflows, includeMaskedBins)) {
0512         const double val = fabs(b.val());
0513         const double vol = includeBinVol? b.dVol() : 1.0;
0514         if (std::isfinite(vol))  ret += val*vol; // aggregate bin volume
0515       }
0516       return ret;
0517     }
0518
0519     /// @brief Convenient alias for areaUnderCurve()
0520     double auc(const bool includeBinVol=true,
0521                const bool includeOverflows=false,
0522                const bool includeMaskedBins=false) const {
0523       return areaUnderCurve(includeBinVol, includeOverflows, includeMaskedBins);
0524     }
0525
0526
0527
0528     /// @brief Construct a covariance matrix from the error breakdown
0529     std::vector<std::vector<double> > covarianceMatrix(const bool ignoreOffDiagonalTerms=false,
0530                                                        const bool includeOverflows=false,
0531                                                        const bool includeMaskedBins=false,
0532                                                        const std::string& pat_uncorr="^stat|^uncor") const {
0533       const size_t nBins = BaseT::numBins(includeOverflows,includeMaskedBins);
0534       std::vector<std::vector<double> > covM(nBins);
0535
0536       // initialise cov matrix to be the right shape
0537       for (size_t i = 0; i < nBins; ++i) {
0538         covM[i] = std::vector<double>(nBins, 0.0);
0539       }
0540
0541       const std::vector<std::string> error_sources = sources();
0542
0543       // nominal-only case, i.e. total uncertainty, labelled as empty string
0544       if (error_sources.size() == 1 && error_sources[0] == "") {
0545         size_t i = 0;
0546         for (const auto& b : BaseT::bins(includeOverflows,includeMaskedBins)) {
0547           covM[i][i] = b.hasSource("")? sqr(0.5*(b.err("").first+b.err("").second)) : 1.0;
0548           ++i;
0549         }
0550         return covM;
0551       }
0552
0553       // more interesting case where we actually have some uncertainty breakdown!
0554       std::smatch match;
0555       const std::regex re(pat_uncorr, std::regex_constants::icase);
0556       for (const std::string& sname : error_sources) {
0557         if (sname == "")  continue;
0558         std::vector<double> systErrs(nBins, 0.0);
0559         size_t k = 0;
0560         for (const auto& b : BaseT::bins(includeOverflows,includeMaskedBins)) {
0561           if (b.hasSource(sname)) {
0562             const auto& errs = b.err(sname); // dn-up pair
0563             systErrs[k] = 0.5 *( fabs(errs.first)+fabs(errs.second));
0564           }
0565           ++k;
0566         }
0567         const bool skipOffDiag = (ignoreOffDiagonalTerms
0568                                   || std::regex_search(sname, match, re));
0569         for (size_t i = 0; i < nBins; ++i) {
0570           for (size_t j = 0; j < nBins; ++j) {
0571             if (skipOffDiag && i != j)  continue;
0572             covM[i][j] += systErrs[i]*systErrs[j];
0573           }
0574         }
0575       }
0576
0577       return covM;
0578     }
0579
0580     /// @}
0581
0582     /// @name Type reductions
0583     /// @{
0584
0585     /// @brief Produce a ScatterND from a EstimateStorage
0586     auto mkScatter(const std::string& path="", const std::string& pat_match = "",
0587                                                const bool includeOverflows=false,
0588                                                const bool includeMaskedBins=false) const {
0589
0590       constexpr size_t N = sizeof...(AxisT);
0591
0592       ScatterND<N+1> rtn;
0593       for (const std::string& a : annotations()) {
0594         if (a != "Type")  rtn.setAnnotation(a, annotation(a));
0595       }
0596       rtn.setAnnotation("Path", path);
0597
0598       const bool incOF = all_CAxes<AxisT...>::value && includeOverflows;
0599       for (const auto& b : BaseT::bins(incOF, includeMaskedBins)) {
0600
0601         // Create the vector of coordinates
0602         Utils::ndarray<double, N+1> vals;
0603         // first fill bin centres, use bin index if axis non-arithmetic
0604         auto indexIfDiscrete = [&vals, &b](auto I) {
0605           using isContinuous = typename BinningT::template is_CAxis<I>;
0606           using isArithmetic = typename BinningT::template is_Arithmetic<I>;
0607           vals[I] = coordPicker<I>(b, std::integral_constant<bool, isArithmetic::value>{},
0608                                       std::integral_constant<bool, isContinuous::value>{});
0609         };
0610         MetaUtils::staticFor<BinningT::Dimension::value>(indexIfDiscrete);
0611         // then fill bin content
0612         vals[N] = b.val();
0613
0614         // Create the vector of error pairs, use 0 if axis not continuous
0615         Utils::ndarray<std::pair<double,double>, N+1> errs;
0616         auto nullifyDiscrete = [&errs, &vals, &b](auto I) {
0617           using isContinuous = typename BinningT::template is_CAxis<I>;
0618           errs[I] = nullifyIfDisc<I>(b, vals[I], std::integral_constant<bool, isContinuous::value>{});
0619         };
0620         MetaUtils::staticFor<BinningT::Dimension::value>(nullifyDiscrete);
0621         const double tot = fabs(b.totalErrPos(pat_match)); // use positive error component
0622         errs[N] = { tot, tot };
0623
0624         // Add the PointND
0625         rtn.addPoint( PointND<N+1>(vals, errs) );
0626       }
0627
0628       // Decorate output scatter with the discrete edges
0629       const BinningT& binning = BaseT::_binning;
0630       auto decorateEdges = [&rtn, &binning](auto I) {
0631         using isContinuous = typename BinningT::template is_CAxis<I>;
0632         if constexpr( !isContinuous::value ) {
0633           const auto& axis  = binning.template axis<I>();
0634           if (axis.numBins()) {
0635             std::stringstream ss;
0636             axis._renderYODA(ss);
0637             rtn.setAnnotation("EdgesA" + std::to_string(I+1), ss.str());
0638           }
0639         }
0640       };
0641       MetaUtils::staticFor<BinningT::Dimension::value>(decorateEdges);
0642
0643       return rtn;
0644     }
0645
0646     /// @brief Split into vector of BinnedEstimates along @a axisN
0647     ///
0648     /// The binning dimension of the returned objects are reduced by one unit.
0649     /// @note Requires at least two binning dimensions.
0650     template<size_t axisN, typename = std::enable_if_t< (axisN < sizeof...(AxisT) && sizeof...(AxisT)>=2) >>
0651     auto mkEstimates(const std::string& path="", const bool includeOverflows=false) const {
0652
0653       // Need to provide a prescription for how to add the two bin contents
0654       auto how2add = [](auto& pivot, const BinType& toCopy) { pivot = toCopy; };
0655       auto rtn = BaseT::template mkBinnedSlices<axisN, BinnedEstimate>(how2add, includeOverflows);
0656       for (const std::string& a : annotations()) {
0657         if (a == "Type")  continue;
0658         for (size_t i = 0; i < rtn.size(); ++i) {
0659           rtn[i].setAnnotation(a, annotation(a));
0660         }
0661       }
0662       for (size_t i = 0; i < rtn.size(); ++i) {
0663         rtn[i].setAnnotation("Path", path);
0664       }
0665
0666       return rtn;
0667     }
0668
0669
0670     /// @brief Method returns clone of the estimate with streamlined error source
0671     AnalysisObject* mkInert(const std::string& path = "",
0672                             const std::string& source = "") const noexcept {
0673       EstimateStorage* rtn = newclone();
0674       rtn->setPath(path);
0675       for (auto& b : rtn->bins(true, true)) {
0676         if (b.numErrs() == 1) {
0677           try {
0678             b.renameSource("", source);
0679           }
0680           catch (YODA::UserError& e) { }
0681         }
0682       }
0683       return rtn;
0684     }
0685
0686     /// @}
0687
0688     /// @name MPI (de-)serialisation
0689     //@{
0690
0691     size_t lengthContent(bool fixed_length = false) const noexcept {
0692       size_t rtn = 0;
0693       for (const auto& bin : BaseT::bins(true, true)) {
0694         rtn += bin._lengthContent(fixed_length);
0695       }
0696       return rtn;
0697     }
0698
0699     std::vector<double> serializeContent(bool fixed_length = false) const noexcept {
0700       std::vector<double> rtn;
0701       const size_t nBins = BaseT::numBins(true, true);
0702       rtn.reserve(nBins * 4);
0703       for (size_t i = 0; i < nBins; ++i) {
0704         const auto& b = BaseT::bin(i);
0705         std::vector<double> bdata = b._serializeContent(fixed_length);
0706         rtn.insert(std::end(rtn),
0707                    std::make_move_iterator(std::begin(bdata)),
0708                    std::make_move_iterator(std::end(bdata)));
0709       }
0710       return rtn;
0711     }
0712
0713
0714     void deserializeContent(const std::vector<double>& data) {
0715
0716       const size_t nBins = BaseT::numBins(true, true);
0717       const size_t minLen = 2*nBins;
0718       if (data.size() < minLen)
0719         throw UserError("Length of serialized data should be at least " + std::to_string(minLen)+"!");
0720
0721       size_t i = 0;
0722       auto itr = data.cbegin();
0723       const auto itrEnd = data.cend();
0724       const bool fixedLen = data.size() == 2*minLen;
0725       while (itr != itrEnd) {
0726         // for estimates, the first element represents the central,
0727         // the subsequent value represents the number of error pairs
0728         const size_t nErrs = fixedLen? 1 : (*(itr + 1) + 0.5); // add 0.5 to avoid rounding issues
0729         auto last = itr + 2*(nErrs+1); // last element + 1
0730         BaseT::bin(i)._deserializeContent(std::vector<double>{itr, last}, fixedLen);
0731         // update for next iteration
0732         itr = last;
0733         ++i;
0734       }
0735     }
0736
0737     // @}
0738
0739   private:
0740
0741     /// @brief Helper function to create a BinningT from
0742     /// a given set @a nBins within a range @a limitsLowUp
0743     template<size_t... Is>
0744     BinningT _mkBinning(const std::vector<size_t>& nBins,
0745                         const std::vector<std::pair<double, double>>& limitsLowUp,
0746                         std::index_sequence<Is...>) const {
0747       return BinningT({((void)Is, Axis<AxisT>(nBins[Is], limitsLowUp[Is].first, limitsLowUp[Is].second))...});
0748     }
0749
0750     /// @brief Helper function to create a BinningT from a scatter @a s
0751     template<size_t... Is>
0752     BinningT _mkBinning(const ScatterND<sizeof...(AxisT)+1>& s, std::index_sequence<Is...>) const {
0753       return BinningT({((void)Is, Axis<AxisT>(s.edges(Is)))...});
0754     }
0755
0756   };
0757
0758
0759
0760   /// @brief User-facing BinnedEstimate class in arbitrary dimensions
0761   template <typename... AxisT>
0762   class BinnedEstimate : public EstimateStorage<AxisT...> {
0763   public:
0764     using EstimateT = BinnedEstimate<AxisT...>;
0765     using BaseT = EstimateStorage<AxisT...>;
0766     using BinType = typename BaseT::BinT;
0767     using Ptr = std::shared_ptr<EstimateT>;
0768
0769     /// @brief Inherit constructors.
0770     using BaseT::BaseT;
0771
0772     BinnedEstimate(const EstimateT&) = default;
0773     BinnedEstimate(EstimateT&&) = default;
0774     BinnedEstimate& operator =(const EstimateT&) = default;
0775     BinnedEstimate& operator =(EstimateT&&) = default;
0776     using AnalysisObject::operator =;
0777
0778     /// @brief Copy constructor (needed for clone functions).
0779     ///
0780     /// @note Compiler won't generate this constructor automatically.
0781     BinnedEstimate(const BaseT& other) : BaseT(other) {}
0782     //
0783     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
0784
0785     /// @brief Move constructor.
0786     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
0787     //
0788     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
0789
0790   };
0791
0792
0793   /// @brief Specialisation of the BinnedEstimate for a 1D histogram
0794   template <typename AxisT>
0795   class BinnedEstimate<AxisT>
0796       : public EstimateStorage<AxisT>,
0797         public XAxisMixin<BinnedEstimate<AxisT>, AxisT> {
0798   public:
0799     using EstimateT = BinnedEstimate<AxisT>;
0800     using BaseT = EstimateStorage<AxisT>;
0801     using BinType = typename BaseT::BinT;
0802     using Ptr = std::shared_ptr<EstimateT>;
0803
0804     /// @brief Inherit constructors.
0805     using BaseT::BaseT;
0806
0807     BinnedEstimate(const EstimateT&) = default;
0808     BinnedEstimate(EstimateT&&) = default;
0809     BinnedEstimate& operator =(const EstimateT&) = default;
0810     BinnedEstimate& operator =(EstimateT&&) = default;
0811     using AnalysisObject::operator =;
0812
0813     /// @brief Copy constructor (needed for clone functions).
0814     ///
0815     /// @note Compiler won't generate this constructor automatically.
0816     BinnedEstimate(const BaseT& other) : BaseT(other) {}
0817     //
0818     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
0819
0820     /// @brief Move constructor.
0821     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
0822     //
0823     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
0824
0825
0826     BinnedEstimate(std::vector<AxisT>&& edges, const std::string& path="", const std::string& title="")
0827               : BaseT(std::move(edges), path, title) {}
0828
0829     BinnedEstimate(const std::vector<AxisT>& edges, const std::string& path="", const std::string& title="")
0830               : BaseT(edges, path, title) {}
0831
0832     /// @brief Constructor with auto-setup of evenly spaced axes.
0833     ///
0834     /// The constructor argument uses double rather than EdgeT to
0835     /// allow for auto-conversion of int to double.
0836     ///
0837     /// @note This constructor is only supported when all axes are continuous.
0838     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT>>
0839     BinnedEstimate(size_t nbins, double lower, double upper, const std::string& path = "", const std::string& title = "")
0840               : BaseT({nbins}, {{lower, upper}}, path, title) {}
0841
0842     /// @brief Make a copy on the stack
0843     EstimateT clone() const noexcept {
0844       return EstimateT(*this);
0845     }
0846
0847     /// @brief Make a copy on the heap
0848     EstimateT* newclone() const noexcept {
0849       return new EstimateT(*this);
0850     }
0851
0852     /// @brief Find bin index for given coordinates
0853     size_t indexAt(const AxisT xCoord) const noexcept {
0854       return BaseT::binAt( {xCoord} ).index();
0855     }
0856
0857     /// @brief Mask/Unmask bin at given set of coordinates
0858     void maskBinAt(const AxisT xCoord, const bool status = true) noexcept {
0859       return BaseT::maskBin({xCoord}, status);
0860     }
0861
0862   };
0863
0864
0865
0866   /// @brief Specialisation of the BinnedEstimate for a 2D BinnedEstimate
0867   template <typename AxisT1, typename AxisT2>
0868   class BinnedEstimate<AxisT1, AxisT2>
0869       : public EstimateStorage<AxisT1, AxisT2>,
0870         public XAxisMixin<BinnedEstimate<AxisT1, AxisT2>, AxisT1>,
0871         public YAxisMixin<BinnedEstimate<AxisT1, AxisT2>, AxisT2> {
0872   public:
0873     using EstimateT = BinnedEstimate<AxisT1, AxisT2>;
0874     using BaseT = EstimateStorage<AxisT1, AxisT2>;
0875     using BinType = typename BaseT::BinT;
0876     using Ptr = std::shared_ptr<EstimateT>;
0877
0878     /// @brief Inherit constructors.
0879     using BaseT::BaseT;
0880
0881     BinnedEstimate(const EstimateT&) = default;
0882     BinnedEstimate(EstimateT&&) = default;
0883     BinnedEstimate& operator =(const EstimateT&) = default;
0884     BinnedEstimate& operator =(EstimateT&&) = default;
0885     using AnalysisObject::operator =;
0886
0887     /// @brief Copy constructor (needed for clone functions).
0888     ///
0889     /// @note Compiler won't generate this constructor automatically.
0890     BinnedEstimate(const BaseT& other) : BaseT(std::move(other)) {}
0891     //
0892     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
0893
0894     /// @brief Move constructor.
0895     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
0896     //
0897     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
0898
0899     BinnedEstimate(std::vector<AxisT1>&& xEdges, std::vector<AxisT2>&& yEdges,
0900                    const std::string& path="", const std::string& title="")
0901               : BaseT(std::move(xEdges), std::move(yEdges), path, title) {}
0902
0903     BinnedEstimate(const std::vector<AxisT1>& xEdges, const std::vector<AxisT2>& yEdges,
0904                    const std::string& path="", const std::string& title="")
0905               : BaseT(xEdges, yEdges, path, title) {}
0906
0907     /// @brief Constructor with auto-setup of evenly spaced axes.
0908     ///
0909     /// The constructor argument uses double rather than EdgeT to
0910     /// allow for auto-conversion of int to double.
0911     ///
0912     /// @note This constructor is only supported when all axes are continuous.
0913     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT1, AxisT2>>
0914     BinnedEstimate(size_t nbinsX, double lowerX, double upperX,
0915                    size_t nbinsY, double lowerY, double upperY,
0916                    const std::string& path = "", const std::string& title = "")
0917               : BaseT({nbinsX, nbinsY}, {{lowerX, upperX}, {lowerY, upperY}}, path, title) {}
0918
0919     /// @brief Make a copy on the stack
0920     EstimateT clone() const noexcept {
0921       return EstimateT(*this);
0922     }
0923
0924     /// @brief Make a copy on the heap
0925     EstimateT* newclone() const noexcept {
0926       return new EstimateT(*this);
0927     }
0928
0929     /// @brief Bin access using global index
0930     BinType& bin(const size_t index) noexcept {
0931       return BaseT::bin(index);
0932     }
0933
0934     /// @brief Bin access using global index (const version)
0935     const BinType& bin(const size_t index) const noexcept {
0936       return BaseT::bin(index);
0937     }
0938
0939     /// @brief Bin access using local indices
0940     BinType& bin(const size_t localX, const size_t localY) noexcept {
0941       return BaseT::bin( {localX, localY} );
0942     }
0943
0944     /// @brief Bin access using local indices (const version)
0945     const BinType& bin(const size_t localX, const size_t localY) const noexcept {
0946       return BaseT::bin( {localX, localY} );
0947     }
0948
0949     /// @brief Bin access using coordinates
0950     BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord) noexcept {
0951       return BaseT::binAt( {xCoord, yCoord} );
0952     }
0953
0954     /// @brief Bin access using coordinates (const version)
0955     const BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord) const noexcept {
0956       return BaseT::binAt( {xCoord, yCoord} );
0957     }
0958
0959     /// @brief Find bin index for given coordinates
0960     size_t indexAt(const AxisT1 xCoord, const AxisT2 yCoord) const noexcept {
0961       return BaseT::binAt( {xCoord, yCoord} ).index();
0962     }
0963
0964     /// @brief Mask/Unmask bin at given set of coordinates
0965     void maskBinAt(const AxisT1 xCoord, const AxisT2 yCoord, const bool status = true) noexcept {
0966       return BaseT::maskBin({xCoord, yCoord}, status);
0967     }
0968
0969   };
0970
0971
0972   /// @brief Specialisation of the BinnedEstimate for a 3D BinnedEstimate
0973   template <typename AxisT1, typename AxisT2, typename AxisT3>
0974   class BinnedEstimate<AxisT1, AxisT2, AxisT3>
0975       : public EstimateStorage<AxisT1, AxisT2, AxisT3>,
0976         public XAxisMixin<BinnedEstimate<AxisT1, AxisT2, AxisT3>, AxisT1>,
0977         public YAxisMixin<BinnedEstimate<AxisT1, AxisT2, AxisT3>, AxisT2>,
0978         public ZAxisMixin<BinnedEstimate<AxisT1, AxisT2, AxisT3>, AxisT3> {
0979   public:
0980     using EstimateT = BinnedEstimate<AxisT1, AxisT2, AxisT3>;
0981     using BaseT = EstimateStorage<AxisT1, AxisT2, AxisT3>;
0982     using BinType = typename BaseT::BinT;
0983     using Ptr = std::shared_ptr<EstimateT>;
0984
0985     /// @brief Inherit constructors.
0986     using BaseT::BaseT;
0987
0988     BinnedEstimate(const EstimateT&) = default;
0989     BinnedEstimate(EstimateT&&) = default;
0990     BinnedEstimate& operator =(const EstimateT&) = default;
0991     BinnedEstimate& operator =(EstimateT&&) = default;
0992     using AnalysisObject::operator =;
0993
0994     /// @brief Copy constructor (needed for clone functions).
0995     ///
0996     /// @note Compiler won't generate this constructor automatically.
0997     BinnedEstimate(const BaseT& other) : BaseT(other) {}
0998     //
0999     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
1000
1001     /// @brief Move constructor.
1002     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
1003     //
1004     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
1005
1006     /// @brief Constructor with auto-setup of evenly spaced axes.
1007     ///
1008     /// The constructor argument uses double rather than EdgeT to
1009     /// allow for auto-conversion of int to double.
1010     ///
1011     /// @note This constructor is only supported when all axes are continuous.
1012     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT1, AxisT2, AxisT3>>
1013     BinnedEstimate(size_t nbinsX, double lowerX, double upperX,
1014                    size_t nbinsY, double lowerY, double upperY,
1015                    size_t nbinsZ, double lowerZ, double upperZ,
1016                    const std::string& path = "", const std::string& title = "")
1017               : BaseT({nbinsX, nbinsY, nbinsZ}, {{lowerX, upperX}, {lowerY, upperY},
1018                       {lowerZ, upperZ}}, path, title) {}
1019
1020     /// @brief Make a copy on the stack
1021     EstimateT clone() const noexcept {
1022       return EstimateT(*this);
1023     }
1024
1025     /// @brief Make a copy on the heap
1026     EstimateT* newclone() const noexcept {
1027       return new EstimateT(*this);
1028     }
1029
1030     /// @brief Bin access using global index
1031     BinType& bin(const size_t index) noexcept {
1032       return BaseT::bin(index);
1033     }
1034
1035     /// @brief Bin access using global index (const version)
1036     const BinType& bin(const size_t index) const noexcept {
1037       return BaseT::bin(index);
1038     }
1039
1040     /// @brief Bin access using local indices
1041     BinType& bin(const size_t localX, const size_t localY, const size_t localZ) noexcept {
1042       return BaseT::bin( {localX, localY, localZ} );
1043     }
1044
1045     /// @brief Bin access using local indices
1046     const BinType& bin(const size_t localX, const size_t localY, const size_t localZ) const noexcept {
1047       return BaseT::bin( {localX, localY, localZ} );
1048     }
1049
1050     /// @brief Bin access using coordinates
1051     BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord) noexcept {
1052       return BaseT::binAt( {xCoord, yCoord, zCoord} );
1053     }
1054
1055     /// @brief Bin access using coordinates (const version)
1056     const BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord) const noexcept {
1057       return BaseT::binAt( {xCoord, yCoord, zCoord} );
1058     }
1059
1060     /// @brief Find bin index for given coordinates
1061     size_t indexAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord) const noexcept {
1062       return BaseT::binAt( {xCoord, yCoord, zCoord} ).index();
1063     }
1064
1065     /// @brief Mask/Unmask bin at given set of coordinates
1066     void maskBinAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord, const bool status = true) noexcept {
1067       return BaseT::maskBin({xCoord, yCoord, zCoord}, status);
1068     }
1069
1070   };
1071
1072   /// @name Combining BinnedEstimates: global operators
1073   /// @{
1074
1075   /// @brief Add two BinnedEstimates
1076   template <typename... AxisT>
1077   inline BinnedEstimate<AxisT...>
1078   operator + (BinnedEstimate<AxisT...> first, const BinnedEstimate<AxisT...>& second) {
1079     first += second;
1080     return first;
1081   }
1082
1083
1084   /// @brief Add two BinnedEstimates
1085   template <typename... AxisT>
1086   inline BinnedEstimate<AxisT...>
1087   operator + (BinnedEstimate<AxisT...> first, BinnedEstimate<AxisT...>&& second) {
1088     first += std::move(second);
1089     return first;
1090   }
1091
1092
1093   /// @brief Subtract two BinnedEstimates
1094   template <typename... AxisT>
1095   inline BinnedEstimate< AxisT...>
1096   operator - (BinnedEstimate<AxisT...> first, const BinnedEstimate<AxisT...>& second) {
1097     first -= second;
1098     return first;
1099   }
1100
1101   /// @brief Subtract two BinnedEstimates
1102   template <typename... AxisT>
1103   inline BinnedEstimate< AxisT...>
1104   operator - (BinnedEstimate<AxisT...> first, BinnedEstimate<AxisT...>&& second) {
1105     first -= std::move(second);
1106     return first;
1107   }
1108
1109   /// @brief Divide two BinnedEstimates
1110   template <typename... AxisT>
1111   inline BinnedEstimate<AxisT...>
1112   divide(const BinnedEstimate<AxisT...>& numer, const BinnedEstimate<AxisT...>& denom,
1113          const std::string& pat_uncorr="^stat|^uncor" ) {
1114     if (numer != denom) {
1115       throw BinningError("Arithmetic operation requires compatible binning!");
1116     }
1117
1118     BinnedEstimate<AxisT...> rtn(numer.binning());
1119     if (numer.path() == denom.path())  rtn.setPath(numer.path());
1120     if (rtn.hasAnnotation("ScaledBy")) rtn.rmAnnotation("ScaledBy");
1121
1122     for (const auto& b_num : numer.bins(true, true)) {
1123       const size_t idx = b_num.index();
1124       rtn.bin(idx) = divide(b_num, denom.bin(idx), pat_uncorr);
1125     }
1126     rtn.maskBins(denom.maskedBins(), true);
1127
1128     return rtn;
1129   }
1130
1131   /// @brief Divide two BinnedEstimates
1132   template <typename... AxisT>
1133   inline BinnedEstimate<AxisT...>
1134   operator / (const BinnedEstimate<AxisT...>& numer, const BinnedEstimate<AxisT...>& denom) {
1135     return divide(numer, denom);
1136   }
1137
1138   /// @brief Divide two BinnedEstimates
1139   template <typename... AxisT>
1140   inline BinnedEstimate<AxisT...>
1141   operator / (BinnedEstimate<AxisT...>&& numer, const BinnedEstimate<AxisT...>& denom) {
1142     return divide(std::move(numer), denom);
1143   }
1144
1145   /// @brief Divide two BinnedEstimates
1146   template <typename... AxisT>
1147   inline BinnedEstimate<AxisT...>
1148   operator / (const BinnedEstimate<AxisT...>& numer, BinnedEstimate<AxisT...>&& denom) {
1149     return divide(numer, std::move(denom));
1150   }
1151
1152   /// @brief Divide two BinnedEstimates
1153   template <typename... AxisT>
1154   inline BinnedEstimate<AxisT...>
1155   operator / (BinnedEstimate<AxisT...>&& numer, BinnedEstimate<AxisT...>&& denom) {
1156     return divide(std::move(numer), std::move(denom));
1157   }
1158
1159
1160   /// @brief Calculate a binned efficiency ratio of two BinnedEstimate objects
1161   ///
1162   /// @note An efficiency is not the same thing as a standard division of two
1163   /// BinnedEstimate objects: the errors are treated as correlated via binomial statistics.
1164   template <typename... AxisT>
1165   inline BinnedEstimate<AxisT...>
1166   efficiency(const BinnedEstimate<AxisT...>& accepted, const BinnedEstimate<AxisT...>& total,
1167              const std::string& pat_uncorr="^stat|^uncor" ) {
1168
1169     if (accepted != total) {
1170       throw BinningError("Arithmetic operation requires compatible binning!");
1171     }
1172
1173     BinnedEstimate<AxisT...> rtn(accepted.binning());
1174
1175     for (const auto& b_acc : accepted.bins(true, true)) {
1176       Estimate est;
1177       const size_t idx = b_acc.index();
1178       try {
1179         est = efficiency(b_acc, total.bin(idx), pat_uncorr);
1180       } catch (const UserError& e) {
1181         //
1182       }
1183       rtn.bin(idx).set(est);
1184     }
1185     return rtn;
1186   }
1187
1188
1189   /// @brief Calculate the asymmetry (a-b)/(a+b) of two BinnedDbn objects
1190   template <typename... AxisT>
1191   inline BinnedEstimate<AxisT...>
1192   asymm(const BinnedEstimate<AxisT...>& a,
1193         const BinnedEstimate<AxisT...>& b,
1194         const std::string& pat_uncorr="^stat|^uncor" ) {
1195     return divde(subtract(a-b, pat_uncorr), add(a+b, pat_uncorr), pat_uncorr);
1196   }
1197
1198
1199   /// @}
1200
1201
1202   /// @name Generalised transformations
1203   /// @{
1204
1205   template <typename... AxisT>
1206   inline void transform(BinnedEstimate<AxisT...>& est, const Trf<1>& fn) {
1207     for (auto& b : est.bins(true, true)) {
1208       b.transform(fn);
1209     }
1210   }
1211
1212   template <typename... AxisT, typename FN>
1213   inline void transform(BinnedEstimate<AxisT...>& est, const FN& fn) {
1214     transform(est, Trf<1>(fn));
1215   }
1216
1217   /// @}
1218
1219
1220   /// @name Convenience aliases
1221   /// @{
1222
1223   /// Define dimension-specific short-hands (Cython sugar)
1224   template<typename A1>
1225   using BinnedEstimate1D = BinnedEstimate<A1>;
1226
1227   template <typename A1, typename A2>
1228   using BinnedEstimate2D = BinnedEstimate<A1, A2>;
1229
1230   template <typename A1, typename A2, typename A3>
1231   using BinnedEstimate3D = BinnedEstimate<A1, A2, A3>;
1232
1233   /// Anonymous namespace to limit visibility
1234   namespace {
1235     template <class T>
1236     struct EstimateMaker;
1237
1238     template<size_t... Is>
1239     struct EstimateMaker<std::index_sequence<Is...>> {
1240       using type = BinnedEstimate< std::decay_t<decltype((void)Is, std::declval<double>())>... >;
1241     };
1242   }
1243
1244   template<size_t N>
1245   using EstimateND = typename EstimateMaker<std::make_index_sequence<N>>::type;
1246
1247
1248   // User-friendly names (continuous axes only)
1249   using Estimate1D = BinnedEstimate<double>;
1250   using Estimate2D = BinnedEstimate<double,double>;
1251   using Estimate3D = BinnedEstimate<double,double,double>;
1252
1253   /// @}
1254
1255 }
1256
1257 #endif