include/YODA/Estimate0D.h

0001 // -*- C++ -*-
0002 //
0003 // This file is part of YODA -- Yet more Objects for Data Analysis
0004 // Copyright (C) 2008-2025 The YODA collaboration (see AUTHORS for details)
0005 //
0006 #ifndef YODA_Estimate0D_h
0007 #define YODA_Estimate0D_h
0008
0009 #include "YODA/AnalysisObject.h"
0010 #include "YODA/Exceptions.h"
0011 #include "YODA/Estimate.h"
0012 #include "YODA/Scatter.h"
0013
0014 #ifdef HAVE_HDF5
0015 #include "YODA/Utils/H5Utils.h"
0016 #endif
0017
0018 #include <string>
0019
0020 namespace YODA {
0021
0022
0023   /// @brief An estimate in 0D
0024   ///
0025   /// This class linkes the Estimate class with
0026   /// AnalysisObject, so it can be used as a type
0027   /// reduction for the Counter type.
0028   class Estimate0D : public Estimate,
0029                      public AnalysisObject {
0030   public:
0031
0032     using Ptr = std::shared_ptr<Estimate0D>;
0033     using BaseT = Estimate;
0034     using BaseT::operator =;
0035     using BaseT::operator +=;
0036     using BaseT::operator -=;
0037     using AnalysisObject::operator =;
0038
0039     /// @name Constructors
0040     /// @{
0041
0042     /// @brief Nullary constructor for unique pointers etc.
0043     ///
0044     /// @note The setting of optional path/title is not possible here in order
0045     /// to avoid overload ambiguity for brace-initialised constructors.
0046     Estimate0D(const std::string& path = "", const std::string& title = "")
0047       : BaseT(), AnalysisObject("Estimate0D", path, title) { }
0048
0049
0050     /// @brief Constructor to set an Estimate0D with a pre-filled state.
0051     ///
0052     /// Principally designed for internal persistency use.
0053     Estimate0D(double v,
0054              std::map<std::string,std::pair<double,double>>& errors,
0055              const std::string& path = "", const std::string& title = "")
0056       : BaseT(v, errors), AnalysisObject("Estimate0D", path, title) { }
0057
0058
0059     /// @brief Alternative constructor to set an Estimate0D with value and uncertainty.
0060     Estimate0D(const double v, const std::pair<double,double>& e, const std::string& source = "",
0061              const std::string& path = "", const std::string& title = "")
0062       : BaseT(v, e, source), AnalysisObject("Estimate0D", path, title) { }
0063
0064     /// @brief Copy constructor (needed for clone functions).
0065     ///
0066     /// @note Compiler won't generate this constructor automatically.
0067     Estimate0D(const Estimate0D& other) : Estimate(other),
0068          AnalysisObject(other.type(), other.path(), other, other.title()) { }
0069
0070     /// @brief Move constructor
0071     Estimate0D(Estimate0D&& other) : BaseT(std::move(other)),
0072          AnalysisObject(other.type(), other.path(), other, other.title()) { }
0073
0074     /// @brief Copy constructor using base class
0075     Estimate0D(const BaseT& other, const std::string& path = "", const std::string& title = "")
0076       : BaseT(other), AnalysisObject("Estimate0D", path, title) { }
0077
0078     /// @brief Move constructor using base class
0079     Estimate0D(BaseT&& other, const std::string& path = "", const std::string& title = "")
0080       : BaseT(std::move(other)), AnalysisObject("Estimate0D", path, title) { }
0081
0082     /// @brief Make a copy on the stack
0083     Estimate0D clone() const noexcept {
0084       return Estimate0D(*this);
0085     }
0086
0087     /// @brief Make a copy on the heap
0088     Estimate0D* newclone() const noexcept {
0089       return new Estimate0D(*this);
0090     }
0091
0092     /// @}
0093
0094
0095     /// @name Operators
0096     /// @{
0097
0098     /// Copy assignment
0099     ///
0100     /// Sets all the parameters using the ones provided from an existing Estimate0D.
0101     Estimate0D& operator=(const Estimate0D& toCopy) noexcept {
0102       if (this != &toCopy) {
0103         AnalysisObject::operator = (toCopy);
0104         BaseT::operator = (toCopy);
0105       }
0106       return *this;
0107     }
0108
0109     /// Move assignment
0110     Estimate0D& operator = (Estimate0D&& toMove) noexcept {
0111       if (this != &toMove) {
0112         AnalysisObject::operator = (toMove);
0113         BaseT::operator = (std::move(toMove));
0114       }
0115       return *this;
0116     }
0117
0118     /// Add two Estimate0Ds
0119     Estimate0D& add(const Estimate0D& toAdd, const std::string& pat_uncorr="^stat|^uncor" ) {
0120       if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
0121       BaseT::add(toAdd, pat_uncorr);
0122       return *this;
0123     }
0124     //
0125     Estimate0D& operator+=(const Estimate0D& toAdd) {
0126       return add(toAdd);
0127     }
0128
0129     /// Add two (rvalue) Estimate0Ds
0130     Estimate0D& add(Estimate0D&& toAdd, const std::string& pat_uncorr="^stat|^uncor" ) {
0131       if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
0132       BaseT::add(std::move(toAdd), pat_uncorr);
0133       return *this;
0134     }
0135     //
0136     Estimate0D& operator+=(Estimate0D&& toAdd) {
0137       return add(std::move(toAdd));
0138     }
0139
0140     /// Subtract two Estimate0Ds
0141     Estimate0D& subtract(const Estimate0D& toSubtract, const std::string& pat_uncorr="^stat|^uncor" ) {
0142       if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
0143       BaseT::subtract(toSubtract, pat_uncorr);
0144       return *this;
0145     }
0146     //
0147     Estimate0D& operator-=(const Estimate0D& toSubtract) {
0148       return subtract(toSubtract);
0149     }
0150
0151     /// Subtract two (rvalue) Estimate0Ds
0152     Estimate0D& subtract(Estimate0D&& toSubtract, const std::string& pat_uncorr="^stat|^uncor" ) {
0153       if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
0154       BaseT::subtract(std::move(toSubtract), pat_uncorr);
0155       return *this;
0156     }
0157     Estimate0D& operator-=(Estimate0D&& toSubtract) {
0158       return subtract(std::move(toSubtract));
0159     }
0160
0161     /// @}
0162
0163     /// @name Dimensions
0164     /// @{
0165
0166     /// @brief Total dimension of this data object
0167     size_t dim() const noexcept { return 1; }
0168
0169     /// @}
0170
0171     /// @name Modifiers
0172     /// @{
0173
0174     /// Reset the internal values
0175     void reset() noexcept {
0176       BaseT::reset();
0177     }
0178
0179     /// @}
0180
0181     /// @name I/O
0182     /// @{
0183
0184     /// @brief Render information about this AO
0185     void _renderYODA(std::ostream& os, const int width = 13) const noexcept {
0186
0187       // Render error sources
0188       const std::vector<std::string> labels = this->sources();
0189       if (labels.size()) {
0190         os << "ErrorLabels: [";
0191         for (size_t i = 0; i < labels.size(); ++i) {
0192           const std::string& src = labels[i];
0193           if (i)  os << ", ";
0194           os << std::quoted(src);
0195         }
0196         os << "]\n";
0197       }
0198
0199       // column header: content types
0200       os << std::setw(width) << std::left << "# value" << "\t";
0201       const int errwidth = std::max(int(std::to_string(labels.size()).size()+7), width); // "errDn(" + src + ")"
0202       for (size_t i = 0; i < labels.size(); ++i) {
0203         const std::string& src = labels[i];
0204         if (src.empty()) {
0205           os << std::setw(errwidth) << std::left << "totalDn" << "\t";
0206           os << std::setw(errwidth) << std::left << "totalUp" << "\t";
0207         }
0208         else {
0209           os << std::setw(errwidth) << std::left << ("errDn(" + std::to_string(i+1) + ")") << "\t";
0210           os << std::setw(errwidth) << std::left << ("errUp(" + std::to_string(i+1) + ")") << "\t";
0211         }
0212       }
0213       os << "\n";
0214
0215       os << std::setw(width) << std::left << val() << "\t"; // render value
0216       // render systs if available
0217       for (const std::string& src : labels) {
0218         if (!hasSource(src)) {
0219           os << std::setw(errwidth) << std::left << "---" << "\t"
0220              << std::setw(errwidth) << std::left << "---" << "\t";
0221           continue;
0222         }
0223         const auto& errs = err(src);
0224         os << std::setw(errwidth) << std::left << errs.first << "\t"
0225            << std::setw(errwidth) << std::left << errs.second << "\t";
0226       }
0227       os << "\n";
0228     }
0229
0230     /// @brief Render scatter-like information about this AO
0231     void _renderFLAT(std::ostream& os, const int width = 13) const noexcept {
0232       const Scatter1D tmp = mkScatter();
0233       tmp._renderYODA(os, width);
0234     }
0235
0236     #ifdef HAVE_HDF5
0237
0238     /// @brief Extract error labels of this AO
0239     void _extractLabels(std::vector<std::string>& labels,
0240                         std::vector<size_t>& labelsizes) const noexcept {
0241       // append new set of keys
0242       std::vector<std::string> keys = sources();
0243       labelsizes.emplace_back(keys.size()+1); //< +1 for length info itself
0244       labels.insert(std::end(labels),
0245                     std::make_move_iterator(std::begin(keys)),
0246                     std::make_move_iterator(std::end(keys)));
0247       // sort and remove duplicates
0248       std::sort(labels.begin(), labels.end());
0249       labels.erase( std::unique(labels.begin(), labels.end()), labels.end() );
0250     }
0251
0252     /// @brief Extract error breakdown of this AO into the map of @a edges
0253     void _extractEdges(std::map<std::string, EdgeHandlerBasePtr>& edges,
0254                        const std::vector<std::string>& labels) const noexcept {
0255
0256       using lenT = EdgeHandler<size_t>;
0257       using lenPtr = EdgeHandlerPtr<size_t>;
0258       const std::string lengthID("sizeinfo");
0259       lenPtr nedges = std::static_pointer_cast<lenT>(edges.find(lengthID)->second);
0260
0261       std::vector<std::string> keys = sources();
0262       nedges->extend({ keys.size() });
0263       const auto& itr = labels.cbegin();
0264       const auto& itrEnd = labels.cend();
0265       for (const string& k : keys) {
0266         size_t dist = std::find(itr, itrEnd, k) - itr;
0267         nedges->extend({ dist });
0268       }
0269
0270     }
0271
0272     #endif
0273
0274     /// @}
0275
0276     /// @name MPI (de-)serialisation
0277     ///@{
0278
0279     size_t lengthContent(bool fixed_length = false) const noexcept {
0280       return BaseT::_lengthContent(fixed_length);
0281     }
0282
0283     std::vector<double> serializeContent(bool fixed_length = false) const noexcept {
0284       return BaseT::_serializeContent(fixed_length);
0285     }
0286
0287     void deserializeContent(const std::vector<double>& data) {
0288       BaseT::_deserializeContent(data, data.size() == 4);
0289     }
0290
0291     /// @}
0292
0293     /// @name Type reductions
0294     /// @{
0295
0296     inline Scatter1D mkScatter(const std::string& path = "",
0297                                const std::string& pat_match = "") const noexcept {
0298       Scatter1D rtn;
0299       for (const std::string& a : annotations()) {
0300         if (a != "Type")  rtn.setAnnotation(a, annotation(a));
0301       }
0302       rtn.setAnnotation("Path", path);
0303
0304       // Add the PointND
0305       const std::pair<double,double> tot = totalErr(pat_match);
0306       rtn.addPoint( Point1D(val(), {fabs(tot.first), tot.second}) );
0307
0308       return rtn;
0309     }
0310
0311     /// @brief Method returns clone of the estimate with streamlined error source
0312     AnalysisObject* mkInert(const std::string& path = "",
0313                             const std::string& source = "") const noexcept {
0314       Estimate0D* rtn = newclone();
0315       rtn->setPath(path);
0316       if (rtn->numErrs() == 1) {
0317         try {
0318           rtn->renameSource("", source);
0319         }
0320         catch (YODA::UserError& e) { }
0321       }
0322       return rtn;
0323     }
0324
0325     /// @}
0326
0327   };
0328
0329   /// @name Generalised transformations
0330   /// @{
0331
0332   inline void transform(Estimate0D& est, const Trf<1>& fn) {
0333     est.transform(fn);
0334   }
0335
0336   template <typename FN>
0337   inline void transform(Estimate0D& est, const FN& fn) {
0338     transform(est, Trf<1>(fn));
0339   }
0340
0341   /// @}
0342
0343
0344   /// @name Global operators for Estimate0D objects
0345   /// @{
0346
0347   /// @brief Add two Estimate0D objects
0348   inline Estimate0D operator + (Estimate0D lhs, const Estimate0D& rhs) {
0349     lhs += rhs;
0350     return lhs;
0351   }
0352
0353   /// @brief Add two Estimate0D objects
0354   inline Estimate0D operator + (Estimate0D lhs, Estimate0D&& rhs) {
0355     lhs += std::move(rhs);
0356     return lhs;
0357   }
0358
0359   /// @brief Subtract two Estimate0D objects
0360   inline Estimate0D operator - (Estimate0D lhs, const Estimate0D& rhs) {
0361     lhs -= rhs;
0362     return lhs;
0363   }
0364
0365   /// @brief Subtract two Estimate0D objects
0366   inline Estimate0D operator - (Estimate0D lhs, Estimate0D&& rhs) {
0367     lhs -= std::move(rhs);
0368     return lhs;
0369   }
0370
0371   /// @brief Divide two Estimate0D objects
0372   inline Estimate0D divide(const Estimate0D& numer, const Estimate0D& denom,
0373                            const std::string& pat_uncorr="^stat|^uncor" ) {
0374     Estimate0D rtn = divide(static_cast<const Estimate&>(numer),
0375                             static_cast<const Estimate&>(denom), pat_uncorr);
0376     if (rtn.hasAnnotation("ScaledBy")) rtn.rmAnnotation("ScaledBy");
0377     if (numer.path() == denom.path())  rtn.setPath(numer.path());
0378     return rtn;
0379   }
0380
0381   /// @brief Divide two Estimate0D objects
0382   inline Estimate0D operator / (const Estimate0D& numer, const Estimate0D& denom) {
0383     return divide(numer, denom);
0384   }
0385
0386   /// @brief Divide two Estimate0D objects
0387   inline Estimate0D operator / (Estimate0D&& numer, const Estimate0D& denom) {
0388     return divide(std::move(numer), denom);
0389   }
0390
0391   /// @brief Divide two Estimate0D objects
0392   inline Estimate0D operator / (const Estimate0D& numer, Estimate0D&& denom) {
0393     return divide(numer, std::move(denom));
0394   }
0395
0396   /// @brief Divide two Estimate0D objects
0397   inline Estimate0D operator / (Estimate0D&& numer, Estimate0D&& denom) {
0398     return divide(std::move(numer), std::move(denom));
0399   }
0400
0401   /// @brief Divide two Estimate0D objects using binomial statistics
0402   inline Estimate0D efficiency(const Estimate0D& accepted, const Estimate0D& total,
0403                                const std::string& pat_uncorr="^stat|^uncor" ) {
0404     Estimate0D rtn = efficiency(static_cast<const Estimate&>(accepted),
0405                                 static_cast<const Estimate&>(total), pat_uncorr);
0406     if (rtn.hasAnnotation("ScaledBy")) rtn.rmAnnotation("ScaledBy");
0407     if (accepted.path() == total.path())  rtn.setPath(total.path());
0408     return rtn;
0409   }
0410
0411   /// @}
0412
0413 }
0414
0415 #endif