EIC code displayed by LXR master/​include/​YODA/​Utils/​MathUtils.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-05-16 08:31:31

 // -*- C++ -*-
 //
 // This file is part of YODA -- Yet more Objects for Data Analysis
 // Copyright (C) 2008-2025 The YODA collaboration (see AUTHORS for details)
 //
 #ifndef YODA_MathUtils_H
 #define YODA_MathUtils_H
 
 /// @todo Add SFINAE math type stuff (see Rivet) and add inrange() and inrange_closed_closed() etc. aliases cf. MCUtils
 
 #include "YODA/Exceptions.h"
 #include "YODA/Config/BuildConfig.h"
 
 #include <algorithm>
 #include <cassert>
 #include <cfloat>
 #include <climits>
 #include <cmath>
 #include <functional>
 #include <iostream>
 #include <limits>
 #include <map>
 #include <numeric>
 #include <numeric>
 #include <ostream>
 #include <sstream>
 #include <stdexcept>
 #include <string>
 #include <utility>
 #include <vector>
 
 namespace YODA {
 
 
   /// Pre-defined numeric type limits
   /// @deprecated Prefer the standard DBL/INT_MAX
   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
 
   /// A pre-defined value of \f$ \pi \f$.
   static const double PI = M_PI;
 
   /// A pre-defined value of \f$ 2\pi \f$.
   static const double TWOPI = 2*M_PI;
 
   /// A pre-defined value of \f$ \pi/2 \f$.
   static const double HALFPI = M_PI_2;
 
   /// Enum for signs of numbers.
   enum Sign { MINUS = -1, ZERO = 0, PLUS = 1 };
 
 
   /// @name Comparison functions for safe (floating point) equality tests
   /// @{
 
   /// @brief Compare a number to zero
   ///
   /// This version for floating point types has a degree of fuzziness expressed
   /// by the absolute @a tolerance parameter, for floating point safety.
   template <typename NUM>
   inline typename std::enable_if_t<std::is_floating_point_v<NUM>, bool>
   isZero(NUM val, double tolerance=1e-5) {
     return fabs(val) < tolerance;
   }
 
   /// @brief Compare a number to zero
   ///
   /// SFINAE template specialisation for integers, since there is no FP
   /// precision issue.
   template <typename NUM>
   inline typename std::enable_if_t<std::is_integral_v<NUM>, bool>
   isZero(NUM val, double = 1e-5) {
     return val==0;
   }
 
   /// @brief Check if a number is NaN
   template <typename NUM>
   inline typename std::enable_if_t<std::is_floating_point_v<NUM>, bool>
   isNaN(NUM val) { return std::isnan(val); }
 
   /// @brief Check if a number is non-NaN
   template <typename NUM>
   inline typename std::enable_if_t<std::is_floating_point_v<NUM>, bool>
   notNaN(NUM val) { return !std::isnan(val); }
 
   /// @brief Compare two numbers for equality with a degree of fuzziness
   ///
   /// This version for floating point types (if any argument is FP) has a degree
   /// of fuzziness expressed by the fractional @a tolerance parameter, for
   /// floating point safety.
   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;
   }
 
   /// @brief Compare two numbers for equality with a degree of fuzziness
   ///
   /// Simpler SFINAE template specialisation for integers, since there is no FP
   /// precision issue.
   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;
   }
 
   /// @brief Comparator wrapper to use with STL algorithms, e.g. std::equal etc.
   static std::function<bool(const double, const double)> fuzzyEqComp =
     [](const double& lhs, const double& rhs) { return fuzzyEquals(lhs, rhs); };
 
 
   /// @brief Compare two numbers for >= with a degree of fuzziness
   ///
   /// The @a tolerance parameter on the equality test is as for @c fuzzyEquals.
   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);
   }
 
   /// @brief Compare two floating point numbers for <= with a degree of fuzziness
   ///
   /// The @a tolerance parameter on the equality test is as for @c fuzzyEquals.
   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);
   }
 
   /// Returns a number floored at the nth decimal place.
   inline double approx(double a, int n = 5) {
     double roundTo = pow(10.0,n);
     a *= roundTo;
     a = floor(a);
     return a/roundTo;
   }
   /// @}
 
 
   /// @name Ranges and intervals
   /// @{
 
   /// Represents whether an interval is open (non-inclusive) or closed (inclusive).
   ///
   /// For example, the interval \f$ [0, \pi) \f$ is closed (an inclusive
   /// boundary) at 0, and open (a non-inclusive boundary) at \f$ \pi \f$.
   enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };
 
 
   /// @brief Determine if @a value is in the range @a low to @a high, for floating point numbers
   ///
   /// Interval boundary types are defined by @a lowbound and @a highbound.
   /// @todo Optimise to one-line at compile time?
   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);
     }
   }
 
   /// Alternative version of inRange for doubles, which accepts a pair for the range arguments.
   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);
   }
 
 
   /// @brief Determine if @a value is in the range @a low to @a high, for integer types
   ///
   /// Interval boundary types are defined by @a lowbound and @a highbound.
   /// @todo Optimise to one-line at compile time?
   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);
     }
   }
 
   /// Alternative version of @c inRange for ints, which accepts a pair for the range arguments.
   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);
   }
 
   /// @}
 
 
   /// @name Miscellaneous numerical helpers
   /// @{
 
   /// Named number-type squaring operation.
   template <typename NUM>
   inline NUM sqr(NUM a) {
     return a*a;
   }
 
   /// Named number-type addition in quadrature operation.
   template <typename... Num>
   inline std::enable_if_t<std::conjunction_v<std::is_arithmetic<Num>...>, double>
   add_quad(Num ... vals) {
     return sqrt((0.0 + ... + sqr(vals)));
   }
 
   /// Find the sign of a number.
   inline int sign(double val) {
     if (isZero(val)) return ZERO;
     const int valsign = (val > 0) ? PLUS : MINUS;
     return valsign;
   }
 
   /// Find the sign of a number.
   inline int sign(int val) {
     if (val == 0) return ZERO;
     return (val > 0) ? PLUS : MINUS;
   }
 
   /// Find the sign of a number.
   inline int sign(long val) {
     if (val == 0) return ZERO;
     return (val > 0) ? PLUS : MINUS;
   }
 
   /// Subtract two numbers with FP fuzziness.
   inline double subtract(double a, double b, double tolerance = 1e-5) {
     if (fuzzyEquals(a,b,tolerance))  return 0.;
     return a - b;
   }
 
   /// Add two numbers with FP fuzziness.
   inline double add(double a, double b, double tolerance = 1e-5) {
     return subtract(a,-b,tolerance);
   }
 
   /// @}
 
 
   /// @name Binning helper functions
   /// @{
 
   /// @brief Make a list of @a nbins + 1 values uniformly spaced between @a xmin and @a xmax inclusive.
   ///
   /// @note The arg ordering and the meaning of the nbins variable is "histogram-like",
   /// as opposed to the Numpy/Matlab version.
   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;
   }
 
 
   /// @brief Make a list of @a nbins + 1 values uniformly spaced in log(x) between @a xmin and @a xmax inclusive.
   ///
   /// @note The arg ordering and the meaning of the nbins variable is "histogram-like",
   /// as opposed to the Numpy/Matlab version, and the xmin and xmax arguments are expressed
   /// in "normal" space, rather than as the logarithms of the xmin/xmax values as in Numpy/Matlab.
   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;
   }
 
 
   /// @todo fspace() for uniform sampling from f(x); requires ability to invert fn... how, in general?
   //inline std::vector<double> fspace(size_t nbins, double xmin, double xmax, std::function<double(double)>& fn)
 
 
   /// @brief Make a list of @a nbins + 1 values spaced with *density* ~ f(x) between @a xmin and @a end inclusive.
   ///
   /// The density function @a fn will be evaluated at @a nsample uniformly
   /// distributed points between @a xmin and @a xmax, its integral approximated
   /// via the Trapezium Rule and used to normalize the distribution, and @a
   /// nbins + 1 edges then selected to (approximately) divide into bins each
   /// containing fraction 1/@a nbins of the integral.
   ///
   /// @note The function @a fn does not need to be a normalized pdf, but it should be non-negative.
   /// Any negative values returned by @a fn(x) will be truncated to zero and contribute nothing to
   /// the estimated integral and binning density.
   ///
   /// @note The arg ordering and the meaning of the nbins variable is "histogram-like",
   /// as opposed to the Numpy/Matlab version, and the xmin and xmax arguments are expressed
   /// in "normal" space, rather than in the function-mapped space as with Numpy/Matlab.
   ///
   /// @note The naming of this function differs from the other, Matlab-inspired ones: the bin spacing is
   /// uniform in the CDF of the density function given, rather than in the function itself. For
   /// most use-cases this is more intuitive.
   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;
   }
 
 
   /// @brief Return the bin index of the given value, @a val, given a vector of bin edges.
   ///
   /// NB. The @a binedges vector must be sorted
   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;
   }
 
   /// @}
 
 }
 
 #endif

This page was automatically generated by the 2.3.7 LXR engine. The LXR team