Estimate.h

Go to the documentation of this file.

// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2023 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_Estimate_h
#define YODA_Estimate_h
 
#include "YODA/Exceptions.h"
#include "YODA/Dbn.h"
#include "YODA/Transformation.h"
#include "YODA/Utils/MathUtils.h"
#include "YODA/Utils/StringUtils.h"
#include <cmath>
#include <map>
#include <string>
#include <regex>
#include <utility>
 
namespace YODA {
 
  class Estimate {
  public:
 
 
 
    Estimate() { reset(); }
 
 
    Estimate(double v, std::map<std::string,std::pair<double,double>>& errors)
      : _value(v), _error(errors) { }
 
 
    Estimate(const double v, const std::pair<double,double>& e, const std::string& source = "")
      : _value(v) { setErr(e, source); }
 
    Estimate(const Estimate& toCopy) {
      _value = toCopy._value;
      _error = toCopy._error;
    }
 
    Estimate(Estimate&& toMove) {
      _value = std::move(toMove._value);
      _error = std::move(toMove._error);
    }
 
    Estimate& operator=(const Estimate& toCopy) noexcept {
      if (this != &toCopy) {
        _value = toCopy._value;
        _error = toCopy._error;
      }
      return *this;
    }
 
    Estimate& operator=(Estimate&& toMove) noexcept {
      if (this != &toMove) {
        _value = std::move(toMove._value);
        _error = std::move(toMove._error);
      }
      return *this;
    }
 
 
 
 
 
    Estimate& add(const Estimate& toAdd, const std::string& pat_uncorr="^stat|^uncor" ) {
      setVal(val() + toAdd.val());
      std::smatch match;
      const std::regex re(pat_uncorr);
      for (const std::string& src : toAdd.sources()) {
        if (!hasSource(src))  setErr(toAdd.errDownUp(src), src);
        else if (std::regex_search(src, match, re)) {
          // treat as uncorrelated between AOs:
          // add in quadrature
          const auto [l_dn,l_up] = errDownUp(src);
          const auto [r_dn,r_up] = toAdd.errDownUp(src);
          const double new_dn = std::sqrt(l_dn*l_dn + r_dn*r_dn);
          const double new_up = std::sqrt(l_up*l_up + r_up*r_up);
          setErr({-new_dn,new_up}, src);
        }
        else {
          // treat as correlated between AOs: add linearly
          const auto [l_dn,l_up] = errDownUp(src);
          const auto [r_dn,r_up] = toAdd.errDownUp(src);
          setErr({l_dn+r_dn, l_up+r_up}, src);
        }
      }
      return *this;
    }
    //
    Estimate& operator+=(const Estimate& toAdd) {
      return add(toAdd);
    }
 
    Estimate& subtract(const Estimate& toSubtract, const std::string& pat_uncorr="^stat|^uncor" ) {
      setVal(val() - toSubtract.val());
      std::smatch match;
      const std::regex re(pat_uncorr);
      for (const std::string& src : toSubtract.sources()) {
        if (!hasSource(src))  setErr(toSubtract.errDownUp(src), src);
        else if (std::regex_search(src, match, re)) {
          // treat as uncorrelated between AOs: add in quadrature
          const auto [l_dn,l_up] = errDownUp(src);
          const auto [r_dn,r_up] = toSubtract.errDownUp(src);
          const double new_dn = std::sqrt(l_dn*l_dn + r_dn*r_dn);
          const double new_up = std::sqrt(l_up*l_up + r_up*r_up);
          setErr({-new_dn,new_up}, src);;
        }
        else {
          // treat as correlated between AOs: add linearly
          const auto [l_dn,l_up] = errDownUp(src);
          const auto [r_dn,r_up] = toSubtract.errDownUp(src);
          setErr({l_dn+r_dn, l_up+r_up}, src);
        }
      }
      return *this;
    }
    //
    Estimate& operator-=(const Estimate& toSubtract) {
      return subtract(toSubtract);
    }
 
 
 
 
    void setValue(const double value) noexcept { _value = value; }
 
    void setVal(const double val) noexcept { setValue(val); }
 
    void setErr(const std::pair<double, double> &err, const std::string& source = "") {
      const std::string& s = Utils::toUpper(source);
      if (s == "TOTAL")
        throw UserError("Use empty string for the total uncertainty!");
      _error[source] = err;
    }
 
    void setErr(const double err, const std::string& source = "") {
      setErr({-fabs(err),fabs(err)}, source);
    }
 
    void set(const double val, const std::pair<double, double> &err, const std::string source = "") {
      setVal(val);
      setErr(err, source);
    }
 
 
    void set(const double val, const double err, const std::string& source = "") {
      setVal(val);
      setErr(err, source);
    }
 
    void reset() noexcept {
      _value = std::numeric_limits<double>::quiet_NaN();
      _error.clear();
    }
 
 
    void scale(const double scalefactor) noexcept {
      _value *= scalefactor;
      for (auto& item : _error) {
        item.second = { item.second.first * scalefactor, item.second.second * scalefactor };
      }
    }
 
    void scale(const Trf<1>& trf) {
      trf.transform(_value, _error);
    }
 
    void transform(const Trf<1>& trf) {
      scale(trf);
    }
 
    void renameSource(const std::string& old_label, const std::string& new_label) {
      if (!hasSource(old_label)) {
        throw UserError("Error map has no such key: "+old_label);
      }
      auto entry = _error.extract(old_label);
      entry.key() = new_label;
      _error.insert(std::move(entry));
    }
 
 
 
  public:
 
 
 
 
    double val() const noexcept { return _value; }
 
    std::pair<double,double> errDownUp(const std::string& source = "") const {
      const size_t count = _error.count(source);
      if (!count)
        throw RangeError("Error map has no such key: "+source);
      return _error.at(source);
    }
 
    std::pair<double,double> err(const std::string& source = "") const {
      return errDownUp(source);
    }
 
    std::pair<double,double> errNegPos(std::string source = "") const {
      return _downUp2NegPos( errDownUp(source) );
    }
 
    double errNeg(const std::string& source = "") const {
      return errNegPos(source).first;
    }
 
    double errPos(const std::string& source = "") const {
      return errNegPos(source).second;
    }
 
    double errDown(const std::string& source = "") const {
      return errDownUp(source).first;
    }
 
    double errUp(const std::string& source = "") const {
      return errDownUp(source).second;
    }
 
    double errAvg(const std::string& source = "") const {
      return _average(errDownUp(source));
    }
 
    std::pair<double,double> relErrDownUp(const std::string& source = "") const {
      auto [dn,up] = errDownUp(source);
      dn = _value != 0 ? dn/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
      up = _value != 0 ? up/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
      return {dn,up};
    }
 
    std::pair<double,double> relErr(const std::string& source = "") const {
      return relErrDownUp(source);
    }
 
    double relErrDown(const std::string& source = "") const {
      return relErr(source).first;
    }
 
    double relErrUp(const std::string& source = "") const {
      return relErr(source).second;
    }
 
    double relErrAvg(const std::string& source = "") const {
      return _average(relErrDownUp(source));
    }
 
    double valMax(const std::string& source = "") const {
      return val() + errPos(source);
    }
 
    double valMin(const std::string& source = "") const {
      return val() + errNeg(source);
    }
 
    std::pair<double,double> quadSum() const noexcept {
      // interpret as { neg , pos } quad sums
      std::pair<double, double> ret = { 0., 0. };
      for (const auto& item : _error) {
        const auto [dn,up] = _downUp2NegPos(item.second);
        ret.first += dn*dn;
        ret.second += up*up;
      }
      return { - std::sqrt(ret.first), std::sqrt(ret.second) };
    }
 
    double quadSumNeg() const noexcept {
      return quadSum().first;
    }
 
    double quadSumPos() const noexcept {
      return quadSum().second;
    }
 
    double quadSumAvg() const noexcept {
      return _average(quadSum());
    }
 
    std::pair<double,double> totalErr() const noexcept {
      // check if the user specified the total uncertainty
      if (_error.count(""))  return _downUp2NegPos(_error.at(""));
      // otherwise return quadrature sum
      return quadSum();
    }
 
    double totalErrNeg() const noexcept {
      return totalErr().first;
    }
 
    double totalErrPos() const noexcept {
      return totalErr().second;
    }
 
    double totalErrAvg() const {
      return _average(totalErr());
    }
 
    double valueErr() const {
      return fabs(totalErrAvg());
    }
 
    std::pair<double,double> relTotalErr() const noexcept {
      auto [neg,pos] = totalErr();
      neg = _value != 0 ? neg/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
      pos = _value != 0 ? pos/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
      return {neg,pos};
    }
 
    double relTotalErrNeg() const noexcept {
      return relTotalErr().first;
    }
 
    double relTotalErrPos() const noexcept {
      return relTotalErr().second;
    }
 
    double relTotalErrAvg() const noexcept {
      return _average(relTotalErr());
    }
 
    std::vector<std::string> sources() const noexcept {
      std::vector<std::string> keys;
      for (const auto& item : _error)  keys.push_back(item.first);
      return keys;
    }
 
    bool hasSource(const std::string& key) const noexcept {
      return _error.count(key);
    }
 
    size_t numErrs() const noexcept {
      return _error.size();
    }
 
 
 
 
    size_t _lengthContent(bool fixed_length = false) const noexcept {
      const size_t nErrs = fixed_length? 1 : numErrs();
      return 2*(nErrs + 1);
    }
 
    std::vector<double> _serializeContent(bool fixed_length = false) const noexcept {
      std::vector<double> rtn;
      const size_t nErrs = fixed_length? 1 : numErrs();
      rtn.reserve(2*(nErrs + 1));
      rtn.push_back(_value);
      if (fixed_length) {
        rtn.push_back(1.0);
        rtn.push_back(totalErrNeg());
        rtn.push_back(totalErrPos());
        return rtn;
      }
      rtn.push_back((double)_error.size());
      for (const auto& item : _error) {
        rtn.push_back(item.second.first);
        rtn.push_back(item.second.second);
      }
      return rtn;
    }
 
    void _deserializeContent(const std::vector<double>& data, bool fixed_length = false) {
 
      if (data.size() < 2)
        throw UserError("Length of serialized data should be at least 2!");
 
      if (2*(fixed_length? 1 : data[1]) != (data.size() - 2))
        throw UserError("Expected "+std::to_string(data[1])+" error pairs!");
 
      reset();
      size_t idx = 0;
      auto itr = data.cbegin();
      const auto itrEnd = data.cend();
      while (itr != itrEnd) {
        if (!idx) {
          _value = *itr; ++itr; ++itr;
        }
        else {
          std::string name("source" + std::to_string(idx));
          const double dn = *itr; ++itr;
          const double up = *itr; ++itr;
          setErr({dn, up}, name);
        }
        ++idx;
      }
      if (numErrs() == 1)  renameSource("source1", "");
    }
 
    std::vector<std::string> serializeSources() const noexcept {
      return sources();
    }
 
    void deserializeSources(const std::vector<std::string>& data) {
 
      const size_t nErrs = numErrs();
      if (data.size() != nErrs)
        throw UserError("Expected " + std::to_string(nErrs) + " error source labels!");
 
      for (size_t i = 0; i < data.size(); ++i) {
        const std::string name("source" + std::to_string(i+1));
        if (!_error.count(name))
          throw UserError("Key names have already been updated!");
        renameSource(name, data[i]);
      }
 
    }
 
    // @}
 
  private:
 
 
 
    std::pair<double,double> _downUp2NegPos(const std::pair<double,double> &e) const noexcept {
      const auto [dn,up] = e;
      if (dn < 0. && up < 0.) {
        // both negative
        const double env = std::min(dn, up);
        return {env, 0.};
      }
      else if (dn < 0.) {
        // dn negative, up positive
        return {dn,up};
      }
      else if (up < 0.) {
        // up negative, dn positive
        return {up, dn};
      }
      // else: both positive
      const double env = std::max(dn, up);
      return {0., env};
    }
 
    double _average(const std::pair<double,double> &err) const noexcept {
      return 0.5*(fabs(err.first) + fabs(err.second));
    }
 
 
 
 
 
    double _value;
 
    std::map<std::string, std::pair<double,double>> _error;
 
 
  public:
 
 
 
    std::string _toString() const noexcept {
      std::string res ="";
      res += ("value="+ std::to_string(_value));
      for (const auto& item : _error) {
        const std::string name = item.first == ""? "user-supplied total " : "";
        res += (", " + name + "error("+ item.first);
        res += (")={ "+ std::to_string(item.second.first));
        res += (", "+ std::to_string(item.second.second)+" }");
      }
      return res;
    }
 
  };
 
 
 
 
  inline Estimate operator + (Estimate lhs, const Estimate& rhs) {
    lhs += rhs;
    return lhs;
  }
 
  inline Estimate operator + (Estimate lhs, Estimate&& rhs) {
    lhs += std::move(rhs);
    return lhs;
  }
 
  inline Estimate operator - (Estimate lhs, const Estimate& rhs) {
    lhs -= rhs;
    return lhs;
  }
 
  inline Estimate operator - (Estimate lhs, Estimate&& rhs) {
    lhs -= std::move(rhs);
    return lhs;
  }
 
  inline Estimate divide(const Estimate& numer, const Estimate& denom,
                         const std::string& pat_uncorr="^stat|^uncor") {
 
    Estimate rtn;
    if (denom.val())  rtn.setVal(numer.val() / denom.val());
    const double newVal = rtn.val();
 
    // get error sources
    std::vector<std::string> sources = numer.sources();
    std::vector<std::string> tmp = denom.sources();
    sources.insert(std::end(sources),
                   std::make_move_iterator(std::begin(tmp)),
                   std::make_move_iterator(std::end(tmp)));
    sources.erase( std::unique(sources.begin(), sources.end()), sources.end() );
 
    std::smatch match;
    const std::regex re(pat_uncorr);
    for (const std::string& src : sources) {
      if (std::regex_search(src, match, re)) {
        // treat as uncorrelated between AOs:
        // add relative errors in quadrature
        double n_dn = 0.0, n_up = 0.0;
        if (numer.hasSource(src)) {
          n_dn = numer.relErrDown(src);
          n_up = numer.relErrUp(src);
        }
        double d_dn = 0.0, d_up = 0.0;
        if (denom.hasSource(src)) {
          d_dn = denom.relErrDown(src);
          d_up = denom.relErrUp(src);
        }
        const double new_dn = fabs(newVal) * std::sqrt(n_dn*n_dn + d_dn*d_dn);
        const double new_up = fabs(newVal) * std::sqrt(n_up*n_up + d_up*d_up);
        rtn.setErr({-new_dn,new_up}, src);
      }
      else {
        // treat as correlated between AOs:
        // work out correlated ratio: R+dR = (N+dN)/(D+dD)
        double n_dn = numer.val(), n_up = numer.val();
        if (numer.hasSource(src)) {
          n_dn += numer.errDown(src);
          n_up += numer.errUp(src);
        }
        double d_dn = denom.val(), d_up = denom.val();
        if (denom.hasSource(src)) {
          d_dn += denom.errDown(src);
          d_up += denom.errUp(src);
        }
        double new_dn = std::numeric_limits<double>::quiet_NaN();
        double new_up = std::numeric_limits<double>::quiet_NaN();
        if (d_dn)  new_dn = n_dn / d_dn - newVal;
        if (d_up)  new_up = n_up / d_up - newVal;
        rtn.setErr({new_dn, new_up}, src);
      }
    }
 
    return rtn;
  }
 
  inline Estimate operator / (const Estimate& numer, const Estimate& denom) {
    return divide(numer, denom);
  }
 
  inline Estimate operator / (Estimate&& numer, const Estimate& denom) {
    return divide(std::move(numer), denom);
  }
 
  inline Estimate operator / (const Estimate& numer, Estimate&& denom) {
    return divide(numer, std::move(denom));
  }
 
  inline Estimate operator / (Estimate&& numer, Estimate&& denom) {
    return divide(std::move(numer), std::move(denom));
  }
 
  inline Estimate efficiency(const Estimate& accepted, const Estimate& total,
                             const std::string& pat_uncorr="^stat|^uncor") {
 
    const double acc_val = accepted.val();
    const double tot_val = total.val();
    if (acc_val > tot_val)
      throw UserError("Attempt to calculate an efficiency when the numerator is not a subset of the denominator: "
                      + Utils::toStr(acc_val) + " / " + Utils::toStr(tot_val));
 
    Estimate rtn = divide(accepted, total, pat_uncorr);
 
    // get error sources
    std::vector<std::string> sources = accepted.sources();
    std::vector<std::string> tmp = total.sources();
    sources.insert(std::end(sources),
                   std::make_move_iterator(std::begin(tmp)),
                   std::make_move_iterator(std::end(tmp)));
    sources.erase( std::unique(sources.begin(), sources.end()), sources.end() );
 
    // set binomial error for uncorrelated error sources
    std::smatch match;
    const double eff = rtn.val();
    const std::regex re(pat_uncorr);
    for (const std::string& src : sources) {
      double err = std::numeric_limits<double>::quiet_NaN();
      if (!tot_val) {
        rtn.setErr({-err,err}, src);
        continue;
      }
      else if (std::regex_search(src, match, re)) {
        const double acc_err = accepted.relTotalErrAvg();
        const double tot_err = total.totalErrAvg();
        err = sqrt(std::abs( ((1-2*eff)*sqr(acc_err) + sqr(eff)*sqr(tot_err)) / sqr(tot_val) ));
        rtn.setErr({-err,err}, src);
        continue;
      }
    }
 
    return rtn;
  }
 
 
}
 
#endif