MathUtils.h

Go to the documentation of this file.

// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2024 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_MathUtils_H
#define YODA_MathUtils_H
 
 
#include "YODA/Exceptions.h"
#include "YODA/Config/BuildConfig.h"
 
#include <algorithm>
#include <functional>
#include <numeric>
#include <cassert>
#include <cfloat>
#include <climits>
#include <cmath>
#include <functional>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <ostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
namespace YODA {
 
 
  const static double MAXDOUBLE = DBL_MAX; // was std::numeric_limits<double>::max(); -- warns in GCC5
  const static double MAXINT = INT_MAX; // was std::numeric_limits<int>::max(); -- warns in GCC5
 
  static const double PI = M_PI;
 
  static const double TWOPI = 2*M_PI;
 
  static const double HALFPI = M_PI_2;
 
  enum Sign { MINUS = -1, ZERO = 0, PLUS = 1 };
 
 
 
 
  template <typename NUM>
  inline typename std::enable_if_t<std::is_floating_point_v<NUM>, bool>
  isZero(NUM val, double tolerance=1e-8) {
    return fabs(val) < tolerance;
  }
 
  template <typename NUM>
  inline typename std::enable_if_t<std::is_integral_v<NUM>, bool>
  isZero(NUM val, double = 1e-5) {
    return val==0;
  }
 
  template <typename NUM>
  inline typename std::enable_if_t<std::is_floating_point_v<NUM>, bool>
  isNaN(NUM val) { return std::isnan(val); }
 
  template <typename NUM>
  inline typename std::enable_if_t<std::is_floating_point_v<NUM>, bool>
  notNaN(NUM val) { return !std::isnan(val); }
 
  template <typename N1, typename N2>
  inline typename std::enable_if_t<
    std::is_arithmetic_v<N1> && std::is_arithmetic_v<N2> &&
   (std::is_floating_point_v<N1> || std::is_floating_point_v<N2>), bool>
  fuzzyEquals(N1 a, N2 b, double tolerance=1e-5) {
    const double absavg = (std::abs(a) + std::abs(b))/2.0;
    const double absdiff = std::abs(a - b);
    const bool rtn = (isZero(a) && isZero(b)) || absdiff < tolerance*absavg;
    return rtn;
  }
 
  template <typename N1, typename N2>
  inline typename std::enable_if_t<
    std::is_integral_v<N1> && std::is_integral_v<N2>, bool>
    fuzzyEquals(N1 a, N2 b, double) { //< NB. unused tolerance parameter for ints, still needs a default value!
    return a == b;
  }
 
  static std::function<bool(const double, const double)> fuzzyEqComp =
    [](const double& lhs, const double& rhs) { return fuzzyEquals(lhs, rhs); };
 
 
  template <typename N1, typename N2>
  inline typename std::enable_if_t<
    std::is_arithmetic_v<N1> && std::is_arithmetic_v<N2>, bool>
  fuzzyGtrEquals(N1 a, N2 b, double tolerance=1e-5) {
    return a > b || fuzzyEquals(a, b, tolerance);
  }
 
  template <typename N1, typename N2>
  inline typename std::enable_if_t<
    std::is_arithmetic_v<N1> && std::is_arithmetic_v<N2>, bool>
  fuzzyLessEquals(N1 a, N2 b, double tolerance=1e-5) {
    return a < b || fuzzyEquals(a, b, tolerance);
  }
 
  inline double approx(double a, int n = 5) {
    double roundTo = pow(10.0,n);
    a *= roundTo;
    a = floor(a);
    return a/roundTo;
  }
 
 
 
 
  enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };
 
 
  template<typename NUM>
  inline bool inRange(NUM value, NUM low, NUM high,
                      RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
    if (lowbound == OPEN && highbound == OPEN) {
      return (value > low && value < high);
    } else if (lowbound == OPEN && highbound == CLOSED) {
      return (value > low && value <= high);
    } else if (lowbound == CLOSED && highbound == OPEN) {
      return (value >= low && value < high);
    } else { // if (lowbound == CLOSED && highbound == CLOSED) {
      return (value >= low && value <= high);
    }
  }
 
  template<typename NUM>
  inline bool inRange(NUM value, std::pair<NUM, NUM> lowhigh,
                      RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
    return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
  }
 
 
  inline bool inRange(int value, int low, int high,
                      RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) {
    if (lowbound == OPEN && highbound == OPEN) {
      return (value > low && value < high);
    } else if (lowbound == OPEN && highbound == CLOSED) {
      return (value > low && value <= high);
    } else if (lowbound == CLOSED && highbound == OPEN) {
      return (value >= low && value < high);
    } else { // if (lowbound == CLOSED && highbound == CLOSED) {
      return (value >= low && value <= high);
    }
  }
 
  inline bool inRange(int value, std::pair<int, int> lowhigh,
                      RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
    return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
  }
 
 
 
 
 
  template <typename NUM>
  inline NUM sqr(NUM a) {
    return a*a;
  }
 
  template <typename Num>
  inline Num add_quad(Num a, Num b) {
    return sqrt(a*a + b*b);
  }
 
  template <typename Num>
  inline Num add_quad(Num a, Num b, Num c) {
    return sqrt(a*a + b*b + c*c);
  }
 
  inline int sign(double val) {
    if (isZero(val)) return ZERO;
    const int valsign = (val > 0) ? PLUS : MINUS;
    return valsign;
  }
 
  inline int sign(int val) {
    if (val == 0) return ZERO;
    return (val > 0) ? PLUS : MINUS;
  }
 
  inline int sign(long val) {
    if (val == 0) return ZERO;
    return (val > 0) ? PLUS : MINUS;
  }
 
 
 
 
 
  inline std::vector<double> linspace(size_t nbins, double xmin, double xmax, bool include_end=true) {
    if (xmax < xmin)  throw RangeError("xmax should not be smaller than xmin!");
    if (nbins == 0)   throw RangeError("Requested number of bins is 0!");
    std::vector<double> rtn;
    const double interval = (xmax-xmin)/static_cast<double>(nbins);
    for (size_t i = 0; i < nbins; ++i) {
      rtn.push_back(xmin + i*interval);
    }
    assert(rtn.size() == nbins);
    if (include_end) rtn.push_back(xmax); // exact xmax, not result of n * interval
    return rtn;
  }
 
 
  inline std::vector<double> logspace(size_t nbins, double xmin, double xmax, bool include_end=true) {
    if (xmax < xmin)  throw RangeError("xmax should not be smaller than xmin!");
    if (xmin < 0)     throw RangeError("xmin should not be negative!");
    if (nbins == 0)   throw RangeError("Requested number of bins is 0!");
    const double logxmin = std::log(xmin);
    const double logxmax = std::log(xmax);
    const std::vector<double> logvals = linspace(nbins, logxmin, logxmax);
    assert(logvals.size() == nbins+1);
    std::vector<double> rtn; rtn.reserve(logvals.size());
    rtn.push_back(xmin);
    for (size_t i = 1; i < logvals.size()-1; ++i) {
      rtn.push_back(std::exp(logvals[i]));
    }
    assert(rtn.size() == nbins);
    if (include_end) rtn.push_back(xmax);
    return rtn;
  }
 
 
  //inline std::vector<double> fspace(size_t nbins, double xmin, double xmax, std::function<double(double)>& fn) {
 
 
  inline std::vector<double> pdfspace(size_t nbins, double xmin, double xmax, std::function<double(double)>& fn, size_t nsample=10000) {
    const double dx = (xmax-xmin)/(double)nsample;
    const std::vector<double> xs = linspace(nsample, xmin, xmax);
    std::vector<double> ys(0, nsample);
    auto posfn = [&](double x){return std::max(fn(x), 0.0);};
    std::transform(xs.begin(), xs.end(), ys.begin(), posfn);
    std::vector<double> areas; areas.reserve(nsample);
    double areasum = 0;
    for (size_t i = 0; i < ys.size()-1; ++i) {
      const double area = (ys[i] + ys[i+1])*dx/2.0;
      areas[i] = area;
      areasum += area;
    }
    const double df = areasum/(double)nbins;
    std::vector<double> xedges{xmin}; xedges.reserve(nbins+1);
    double fsum = 0;
    for (size_t i = 0; i < nsample-1; ++i) {
      fsum += areas[i];
      if (fsum > df) {
        fsum = 0;
        xedges.push_back(xs[i+1]);
      }
    }
    xedges.push_back(xmax);
    assert(xedges.size() == nbins+1);
    return xedges;
  }
 
 
  template <typename NUM>
  inline int index_between(const NUM& val, const std::vector<NUM>& binedges) {
    if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range
    int index = -1;
    for (size_t i = 1; i < binedges.size(); ++i) {
      if (val < binedges[i]) {
        index = i-1;
        break;
      }
    }
    assert(inRange(index, -1, binedges.size()-1));
    return index;
  }
 
 
 
 
 
  inline double effNumEntries(const double sumW, const double sumW2) {
    if (isZero(sumW2))  return 0;
    return sqr(sumW) / sumW2;
  }
 
  inline double effNumEntries(const std::vector<double>& weights) {
    double sumW = 0.0, sumW2 = 0.0;
    for (size_t i = 0; i < weights.size(); ++i) {
      sumW += weights[i];
      sumW2 += sqr(weights[i]);
    }
    return effNumEntries(sumW, sumW2);
  }
 
  inline double mean(const std::vector<int>& sample) {
    double mean = 0.0;
    for (size_t i=0; i<sample.size(); ++i) {
      mean += sample[i];
    }
    return mean/sample.size();
  }
 
  inline double mean(const double sumWX, const double sumW) {
    return sumW? sumWX / sumW : std::numeric_limits<double>::quiet_NaN();
  }
 
  inline double mean(const std::vector<double>& sample,
                     const std::vector<double>& weights) {
    if (sample.size() != weights.size())  throw RangeError("Inputs should have equal length!");
    double sumWX = 0., sumW = 0.;
    for (size_t i = 0; i < sample.size(); ++i) {
      sumW  += weights[i];
      sumWX += weights[i]*sample[i];
    }
    return mean(sumWX, sumW);
  }
 
  inline double variance(const double sumWX, const double sumW,
                         const double sumWX2, const double sumW2) {
    const double num = sumWX2*sumW - sqr(sumWX);
    const double den = sqr(sumW) - sumW2;
    // if (fabs(num) < 1e-10 && fabs(den) < 1e-10) {
    //   return std::numeric_limits<double>::quiet_NaN();
    // }
    return den? fabs(num/den): std::numeric_limits<double>::quiet_NaN();
  }
 
  inline double variance(const std::vector<double>& sample,
                         const std::vector<double>& weights) {
    if (sample.size() != weights.size())  throw RangeError("Inputs should have equal length!");
    if (fuzzyLessEquals(effNumEntries(weights), 1.0)) {
       //throw LowStatsError("Requested variance of a distribution with only one effective entry");
       return std::numeric_limits<double>::quiet_NaN();
    }
    double sumWX = 0., sumW = 0.;
    double sumWX2 = 0., sumW2 = 0.;
    for (size_t i = 0; i < sample.size(); ++i) {
      sumW   += weights[i];
      sumWX  += weights[i]*sample[i];
      sumW2  += sqr(weights[i]);
      sumWX2 += weights[i]*sqr(sample[i]);
    }
    return variance(sumWX, sumW, sumWX2, sumW2);
  }
 
  inline double stdDev(const double sumWX, const double sumW,
                       const double sumWX2, const double sumW2) {
    return std::sqrt(variance(sumWX, sumW, sumWX2, sumW2));
  }
 
  inline double stdDev(const std::vector<double>& sample,
                       const std::vector<double>& weights) {
    return std::sqrt(variance(sample, weights));
  }
 
  inline double stdErr(const double sumWX, const double sumW,
                       const double sumWX2, const double sumW2) {
    const double effN = effNumEntries(sumW, sumW2);
    if (effN == 0)  return std::numeric_limits<double>::quiet_NaN();
    const double var = variance(sumWX, sumW, sumWX2, sumW2);
    return std::sqrt(var / effN);
  }
 
  inline double stdErr(const std::vector<double>& sample,
                        const std::vector<double>& weights) {
    if (sample.size() != weights.size())  throw RangeError("Inputs should have equal length!");
    const double effN = effNumEntries(weights);
    if (effN == 0)  return std::numeric_limits<double>::quiet_NaN();
    const double var = variance(sample, weights);
    return std::sqrt(var / effN);
  }
 
  inline double RMS(const double sumWX2, const double sumW, const double sumW2) {
    // Weighted RMS defined as
    // rms = sqrt(sum{w x^2} / sum{w})
    const double effN = effNumEntries(sumW, sumW2);
    if (effN == 0)  return std::numeric_limits<double>::quiet_NaN();
    const double meanSq = sumWX2 / sumW;
    return std::sqrt(meanSq);
  }
 
  inline double RMS(const std::vector<double>& sample,
                    const std::vector<double>& weights) {
    if (sample.size() != weights.size())  throw RangeError("Inputs should have equal length!");
    double sumWX2 = 0., sumW = 0., sumW2 = 0.;
    for (size_t i = 0; i < sample.size(); ++i) {
      sumW   += weights[i];
      sumW2  += sqr(weights[i]);
      sumWX2 += weights[i]*sqr(sample[i]);
    }
    return RMS(sumWX2, sumW, sumW2);
  }
 
  inline double covariance(const std::vector<int>& sample1, const std::vector<int>& sample2) {
    const double mean1 = mean(sample1);
    const double mean2 = mean(sample2);
    const size_t N = sample1.size();
    double cov = 0.0;
    for (size_t i = 0; i < N; i++) {
      const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
      cov += cov_i;
    }
    if (N > 1) return cov/(N-1);
    else return 0.0;
  }
 
 
  inline double correlation(const std::vector<int>& sample1, const std::vector<int>& sample2) {
    const double cov = covariance(sample1, sample2);
    const double var1 = covariance(sample1, sample1);
    const double var2 = covariance(sample2, sample2);
    const double correlation = cov/sqrt(var1*var2);
    const double corr_strength = correlation*sqrt(var2/var1);
    return corr_strength;
  }
 
 
  inline double naiveChi2(const std::vector<double>& sample1, const std::vector<double>& sample2,
                     const std::vector<double>& s1errors = std::vector<double>{},
                     const std::vector<double>& s2errors = std::vector<double>{}) {
    if (sample1.size() != sample2.size()) {
      throw RangeError("Inputs should have equal length!");
    }
    if (s1errors.size() && sample1.size() != s1errors.size()) {
      throw RangeError("Inputs should have equal length!");
    }
    if (s2errors.size() && sample2.size() != s2errors.size()) {
      throw RangeError("Inputs should have equal length!");
    }
    const size_t N = sample1.size();
    double chi2 = 0.0;
    for (size_t i = 0; i < N; ++i) {
      double temp = sqr(sample1[i] - sample2[i]);
      if (s1errors.size()) {
        temp /= sqr(s1errors[i]) + sqr(s2errors[i]);
      }
      chi2 += temp;
    }
    return chi2;
  }
 
  inline double naiveChi2reduced(const std::vector<double>& sample1, const std::vector<double>& sample2,
                        const std::vector<double>& s1errors = std::vector<double>{},
                        const std::vector<double>& s2errors = std::vector<double>{}) {
    if (sample1.empty()) throw RangeError("Inputs should not have 0 length!");
    return naiveChi2(sample1, sample2, s1errors, s2errors)/sample1.size();
  }
 
 
 
}
 
#endif