MathFunctions.h

Go to the documentation of this file.

// Copyright (c) 2002-present, The OpenMS Team -- 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 <boost/random/mersenne_twister.hpp> // for mt19937_64
#include <boost/random/uniform_int.hpp>
#include <cmath>
#include <utility> // for std::pair
 
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;
  }
 
  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)
  {
    if (x >= T(0)) { return T(floor(x + T(0.5))); }
    else { return T(ceil(x - T(0.5))); }
  }
 
  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 ppm Parts-per-million error
    @param 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);
    }
  };
} // namespace Math
} // namespace OpenMS