MathFunctions.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg$
// $Authors: Marc Sturm $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/CONCEPT/Macros.h>
#include <OpenMS/CONCEPT/Types.h>
#include <OpenMS/KERNEL/RangeManager.h>
 
#include <boost/random/mersenne_twister.hpp> // for mt19937_64
#include <boost/random/uniform_int.hpp>
#include <cmath>
#include <boost/math/special_functions/binomial.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/distributions/binomial.hpp>
#include <boost/math/distributions/complement.hpp>
#include <limits>
#include <utility> // for std::pair
#include <vector>
 
namespace OpenMS
{
namespace Math
{
 
  template<typename T>
  bool extendRange(T& min, T& max, const T& value)
  {
    if (value < min)
    {
      min = value;
      return true;
    }
    if (value > max)
    {
      max = value;
      return true;
    }
    return false;
  }
 
  template<typename T>
  bool contains(T value, T min, T max)
  {
    return min <= value && value <= max;
  }
 
  inline std::pair<double, double> zoomIn(const double left, const double right, const float factor, const float align)
  {
    OPENMS_PRECONDITION(factor >= 0, "Factor must be >=0")
    OPENMS_PRECONDITION(align >= 0, "align must be >=0")
    OPENMS_PRECONDITION(align <= 1, "align must be <=1")
    std::pair<double, double> res;
    auto old_width = right - left;
    auto offset_left = (1.0f - factor) * old_width * align;
    res.first = left + offset_left;
    res.second = res.first + old_width * factor;
    return res;
  }
 
  using BinContainer = std::vector<RangeBase>;
  inline BinContainer createBins(double min, double max, uint32_t number_of_bins, double extend_margin = 0)
  {
    OPENMS_PRECONDITION(number_of_bins >= 1, "Number of bins must be >= 1")
    OPENMS_PRECONDITION(min < max, "Require min < max");
    std::vector<RangeBase> res(number_of_bins);
    const double bin_width = (max - min) / number_of_bins;
    for (uint32_t i = 0; i < number_of_bins; ++i)
    {
      res[i] = RangeBase(min + i * bin_width, min + (i + 1) * bin_width);
      res[i].extendLeftRight(extend_margin);
    }
    res.front().setMin(min); // undo potential margin
    res.back().setMax(max);  // undo potential margin
    
    return res;
  }
 
 
  inline double ceilDecimal(double x, int decPow)
  {
    return (ceil(x / pow(10.0, decPow))) * pow(10.0, decPow); // decimal shift right, ceiling, decimal shift left
  }
 
  inline double roundDecimal(double x, int decPow)
  {
    if (x > 0) return (floor(0.5 + x / pow(10.0, decPow))) * pow(10.0, decPow);
 
    return -((floor(0.5 + fabs(x) / pow(10.0, decPow))) * pow(10.0, decPow));
  }
 
  inline double intervalTransformation(double x, double left1, double right1, double left2, double right2)
  {
    return left2 + (x - left1) * (right2 - left2) / (right1 - left1);
  }
 
  inline double linear2log(double x)
  {
    return log10(x + 1); //+1 to avoid negative logarithms
  }
 
  inline double log2linear(double x)
  {
    return pow(10, x) - 1;
  }
 
  inline bool isOdd(UInt x)
  {
    return (x & 1) != 0;
  }
 
  template<typename T>
  T round(T x)
  {
    return std::round(x);
  }
 
  template<typename T>
  T roundTo(const T value, int digits)
  {
    T factor = 1.0;
    if (digits > 0)
    {
      for (int i = 0; i < digits; ++i)
        factor *= 10.0;
    }
    else if (digits < 0)
    {
      for (int i = 0; i < -digits; ++i)
        factor /= 10.0;
    }
 
    return std::round(value * factor) / factor;
  }
 
  template<typename T>
  double percentOf(T value, T total, int digits)
  {
    if (value < 0) { throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Value must be non-negative", String(value)); }
    if (total < 0) { throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Total must be non-negative", String(total)); }
    if (total <= 0) // avoid float equality compare
    {
      return 0.0; // avoid division by zero
    }
    return roundTo(value * 100.0 / total, digits);
  }
 
  inline bool approximatelyEqual(double a, double b, double tol)
  {
    return std::fabs(a - b) <= tol;
  }
 
  template<typename T>
  T gcd(T a, T b)
  {
    T c;
    while (b != 0)
    {
      c = a % b;
      a = b;
      b = c;
    }
    return a;
  }
 
  template<typename T>
  T gcd(T a, T b, T& u1, T& u2)
  {
    u1 = 1;
    u2 = 0;
    T u3 = a;
 
    T v1 = 0;
    T v2 = 1;
    T v3 = b;
 
    while (v3 != 0)
    {
      T q = u3 / v3;
      T t1 = u1 - v1 * q;
      T t2 = u2 - v2 * q;
      T t3 = u3 - v3 * q;
 
      u1 = v1;
      u2 = v2;
      u3 = v3;
 
      v1 = t1;
      v2 = t2;
      v3 = t3;
    }
 
    return u3;
  }
 
  template<typename T>
  T getPPM(T mz_obs, T mz_ref)
  {
    return (mz_obs - mz_ref) / mz_ref * 1e6;
  }
 
  template<typename T>
  T getPPMAbs(T mz_obs, T mz_ref)
  {
    return std::fabs(getPPM(mz_obs, mz_ref));
  }
 
  template<typename T>
  T ppmToMass(T ppm, T mz_ref)
  {
    return (ppm / T(1e6)) * mz_ref;
  }
 
  /*
    @brief Compute the absolute mass diff in [Th], given a ppm value and a reference point.
 
    The returned mass diff is always positive!
 
    @param[in] ppm Parts-per-million error
    @param[in] mz_ref Reference m/z
    @return The absolute mass diff in [Th]
  */
  template<typename T>
  T ppmToMassAbs(T ppm, T mz_ref)
  {
    return std::fabs(ppmToMass(ppm, mz_ref));
  }
 
  inline std::pair<double, double> getTolWindow(double val, double tol, bool ppm)
  {
    double left, right;
 
    if (ppm)
    {
      left = val - val * tol * 1e-6;
      right = val / (1.0 - tol * 1e-6);
    }
    else
    {
      left = val - tol;
      right = val + tol;
    }
 
    return std::make_pair(left, right);
  }
 
  template<typename T1>
  typename T1::value_type quantile(const T1& x, double q)
  {
    if (x.empty()) throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Quantile requested from empty container.");
    if (q < 0.0) q = 0.;
    if (q > 1.0) q = 1.;
 
    const auto n = x.size();
    const auto id = std::max(0., n * q - 1); // -1 for c++ index starting at 0
    const auto lo = floor(id);
    const auto hi = ceil(id);
    const auto qs = x[lo];
    const auto h = (id - lo);
 
    return (1.0 - h) * qs + h * x[hi];
  }
 
  // portable random shuffle
  class OPENMS_DLLAPI RandomShuffler
  {
  public:
    explicit RandomShuffler(int seed): rng_(boost::mt19937_64(seed))
    {
    }
 
    explicit RandomShuffler(const boost::mt19937_64& mt_rng): rng_(mt_rng)
    {
    }
 
    RandomShuffler() = default;
    ~RandomShuffler() = default;
 
    boost::mt19937_64 rng_;
    template<class RandomAccessIterator>
    void portable_random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
    {
      for (auto i = (last - first) - 1; i > 0; --i) // OMS_CODING_TEST_EXCLUDE
      {
        boost::uniform_int<decltype(i)> d(0, i);
        std::swap(first[i], first[d(rng_)]);
      }
    }
 
    void seed(uint64_t val)
    {
      rng_.seed(val);
    }
  };
 
  inline double log_binomial_coef(unsigned n, unsigned k) 
  {
    // Handle edge cases for improved numerical stability
    if (k > n) 
    {
      throw std::invalid_argument("k cannot be greater than n in binomial coefficient");
    }
    
    if (k == 0 || k == n) 
    {
      return 0.0;  // log(1) = 0
    }
    
    // Use symmetry to minimize computation for large k
    if (k > n / 2) 
    {
      k = n - k;
    }
    
    return boost::math::lgamma(n + 1.0) - boost::math::lgamma(k + 1.0) - boost::math::lgamma(n - k + 1.0);
  }
 
  inline double log_sum_exp(double x, double y) 
  {
    // Handle infinite cases
    if (std::isinf(x) && x < 0) return y;
    if (std::isinf(y) && y < 0) return x;
    
    // Use the maximum value for numerical stability
    double max_val = std::max(x, y);
    return max_val + std::log(std::exp(x - max_val) + std::exp(y - max_val));
  }
 
  inline double binomial_cdf_complement(unsigned N, unsigned n, double p)
  {
    if (p < 0.0 || p > 1.0)
    {
      throw std::invalid_argument("Probability p must be between 0 and 1");
    }
    if (n > N)
    {
      throw std::invalid_argument("n cannot be greater than N");
    }
 
    if (n == 0)   return 1.0;                // P(X ≥ 0) = 1
    if (p == 0.0) return (n == 0) ? 1.0 : 0.0;
    if (p == 1.0) return 1.0;               // all mass at N
 
    const boost::math::binomial_distribution<double> dist(N, p);
    return boost::math::cdf(boost::math::complement(dist, n - 1));
  }
} // namespace Math
} // namespace OpenMS