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