Tools/CIC/utl/Math.h

Go to the documentation of this file.

// \author Darko Veberic
// \date 2014

#ifndef _cic_utl_Math_h_
#define _cic_utl_Math_h_

#include <TRandom3.h>
#include <TMath.h>
#include <cmath>
#include <functional>
#include <iostream>
#include <algorithm>
#include <numeric>


namespace utl {

  constexpr double degree = M_PI / 180;


  template<typename T>
  inline
  constexpr
  T
  Sqr(const T& x)
  {
    return x*x;
  }


  inline
  double
  Sin2(const double x)
  {
    return Sqr(std::sin(x));
  }


  inline
  double
  Sin2Deg(const double x)
  {
    return Sin2(x*degree);
  }


  extern TRandom3* const gRandom3;


  inline
  double
  Fermi(const double x)
  {
    return 1 / (1 + std::exp(-x));
  }


  inline
  constexpr
  double
  StepFunction(const double x)
  {
    return x > 0;
  }


  inline
  constexpr
  double
  LinearSigmoid(const double x, const double delta)
  {
    return (x < -delta) ?
      0 :
      ((x > delta) ?
        1 :
        0.5*(1 + x/delta));
  }


  inline
  double
  QuadraticSigmoid(const double x, const double delta)
  {
    // for x in interval [-2*delta, 2*delta] gives smooth transiton from 0 to 1.
    const double d = 2*delta;
    return
      (x <= -d) ?
        0 :
        ((x >= d) ?
           1 :
           ((x <= 0) ?
              0.5*Sqr(1 + x/d) :
              1 - 0.5*Sqr(x/d - 1)
           )
        );
  }


  const std::function<double(double, double)> gQuad =
    [](const double x, const double d) { return QuadraticSigmoid(x, d); };


  inline
  double
  ErfSigmoid(const double x, const double delta)
  {
    return 0.5*(1 + std::erf(x/delta));
  }


  inline
  constexpr
  double
  Interpolate(const double dx, const double dy,
              const double x)
  {
    return dy * x / dx;
  }


  inline
  constexpr
  double
  Interpolate(const double x1, const double y1,
              const double x2, const double y2,
              const double x)
  {
    return y1 + Interpolate(x2 - x1, y2 - y1, x - x1);
  }


  // Kolmogorov-Smirnov measure of uniformity
  double UniformKolmogorovSmirnov(const std::vector<double>& cdf);

  double UniformKolmogorovSmirnov(const std::vector<std::pair<double, double>>& wcdf,
                                  const double totalWeight = 0);


  // Anderson-Darling measure of uniformity
  double UniformAndersonDarling(const std::vector<std::pair<double, double>>& wcdf,
                                const double totalWeight = 0);


  inline
  double
  UniformAndersonDarling(const std::vector<std::vector<std::pair<double, double>>>& wcdf)
  {
    double ad = 0;
    for (const auto& row : wcdf)
      ad += UniformAndersonDarling(row);
    return ad;
  }


  // Chi^2 measure of uniformity
  template<typename T>
  double
  UniformChi2Ndof(const std::vector<T>& counts, double nTot = 0)
  {
    if (!nTot) {
      nTot = std::accumulate(counts.begin(), counts.end(), T(0));
      if (!nTot)
        return 0;
    }
    const unsigned int n = counts.size();
    const double nExp = nTot / n;
    double chi2 = 0;
    for (const T c : counts)
      chi2 += Sqr(nExp - c);
    return chi2 / (nExp * (n - 1));
  }


  template<typename T>
  double
  UniformChi2Ndof(const std::vector<std::vector<T>>& counts)
  {
    double chi2 = 0;
    double nTot = 0;
    for (const auto& row : counts) {
      const double rowTot = std::accumulate(row.begin(), row.end(), T(0));
      if (!rowTot)
        continue;
      nTot += rowTot;
      const unsigned int bins = row.size();
      const double mean = rowTot / bins;
      double rowChi2 = 0;
      for (const T nBin : row)
        rowChi2 += Sqr(mean - nBin);
      chi2 += rowChi2 * bins / (bins - 1);
    }
    return chi2 / nTot;
  }


  inline
  double
  PoissonDeviance(const double mean, const double n)
  {
    using namespace std;
    if (mean < 0 || n < 0)
      return numeric_limits<double>::infinity();
    if (mean == 0)
      return (n == 0) ? 0 : numeric_limits<double>::infinity();
    return mean - n * log(mean) + TMath::LnGamma(n + 1);
  }


  inline
  double
  PoissonProbability(const double mean, const double n)
  {
    using namespace std;
    if (mean < 0 || n < 0)
      return 0;
    if (mean == 0)
      return n == 0;
    return std::exp(-PoissonDeviance(mean, n));
  }


  inline
  double
  NemesLnGamma(const double z)
  {
    return 0.5*(std::log(2*M_PI) - std::log(z)) + z * (std::log(z + 1/(12*z - 1/(10*z))) - 1);
  }


  template<typename T>
  double
  UniformPoissonChi2(const std::vector<T>& counts, double nTot = 0)
  {
    if (!nTot) {
      nTot = std::accumulate(counts.begin(), counts.end(), T(0));
      if (!nTot)
        return 0;
    }
    const unsigned int n = counts.size();
    const double nExp = nTot / n;
    double chi2 = 0;
    for (const T c : counts)
      chi2 += (c ? c * std::log(c) : 0);
    chi2 -= nTot * std::log(nExp);
    return 2*chi2;
  }


  template<typename T>
  inline
  double
  UniformPoissonChi2(const std::vector<std::vector<T>>& counts)
  {
    double chi2 = 0;
    for (const auto& row : counts)
      chi2 += UniformPoissonChi2(row);
    return chi2;
  }

}


#endif