OpenMS
MathFunctions.h
Go to the documentation of this file.
1 // Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
2 // SPDX-License-Identifier: BSD-3-Clause
3 //
4 // --------------------------------------------------------------------------
5 // $Maintainer: Timo Sachsenberg$
6 // $Authors: Marc Sturm $
7 // --------------------------------------------------------------------------
8 
9 #pragma once
10 
12 #include <OpenMS/CONCEPT/Macros.h>
13 #include <OpenMS/CONCEPT/Types.h>
15 
16 #include <boost/random/mersenne_twister.hpp> // for mt19937_64
17 #include <boost/random/uniform_int.hpp>
18 #include <cmath>
19 #include <utility> // for std::pair
20 #include <vector>
21 
22 namespace OpenMS
23 {
31 namespace Math
32 {
33 
42  template<typename T>
43  bool extendRange(T& min, T& max, const T& value)
44  {
45  if (value < min)
46  {
47  min = value;
48  return true;
49  }
50  if (value > max)
51  {
52  max = value;
53  return true;
54  }
55  return false;
56  }
57 
63  template<typename T>
64  bool contains(T value, T min, T max)
65  {
66  return min <= value && value <= max;
67  }
68 
89  inline std::pair<double, double> zoomIn(const double left, const double right, const float factor, const float align)
90  {
91  OPENMS_PRECONDITION(factor >= 0, "Factor must be >=0")
92  OPENMS_PRECONDITION(align >= 0, "align must be >=0")
93  OPENMS_PRECONDITION(align <= 1, "align must be <=1")
94  std::pair<double, double> res;
95  auto old_width = right - left;
96  auto offset_left = (1.0f - factor) * old_width * align;
97  res.first = left + offset_left;
98  res.second = res.first + old_width * factor;
99  return res;
100  }
101 
102  using BinContainer = std::vector<RangeBase>;
118  inline BinContainer createBins(double min, double max, uint32_t number_of_bins, double extend_margin = 0)
119  {
120  OPENMS_PRECONDITION(number_of_bins >= 1, "Number of bins must be >= 1")
121  OPENMS_PRECONDITION(min < max, "Require min < max");
122  std::vector<RangeBase> res(number_of_bins);
123  const double bin_width = (max - min) / number_of_bins;
124  for (uint32_t i = 0; i < number_of_bins; ++i)
125  {
126  res[i] = RangeBase(min + i * bin_width, min + (i + 1) * bin_width);
127  res[i].extendLeftRight(extend_margin);
128  }
129  res.front().setMin(min); // undo potential margin
130  res.back().setMax(max); // undo potential margin
131 
132  return res;
133  }
134 
135 
147  inline double ceilDecimal(double x, int decPow)
148  {
149  return (ceil(x / pow(10.0, decPow))) * pow(10.0, decPow); // decimal shift right, ceiling, decimal shift left
150  }
151 
162  inline double roundDecimal(double x, int decPow)
163  {
164  if (x > 0) return (floor(0.5 + x / pow(10.0, decPow))) * pow(10.0, decPow);
165 
166  return -((floor(0.5 + fabs(x) / pow(10.0, decPow))) * pow(10.0, decPow));
167  }
168 
174  inline double intervalTransformation(double x, double left1, double right1, double left2, double right2)
175  {
176  return left2 + (x - left1) * (right2 - left2) / (right1 - left1);
177  }
178 
186  inline double linear2log(double x)
187  {
188  return log10(x + 1); //+1 to avoid negative logarithms
189  }
190 
198  inline double log2linear(double x)
199  {
200  return pow(10, x) - 1;
201  }
202 
208  inline bool isOdd(UInt x)
209  {
210  return (x & 1) != 0;
211  }
212 
218  template<typename T>
219  T round(T x)
220  {
221  return std::round(x);
222  }
223 
238  template<typename T>
239  T roundTo(const T value, int digits)
240  {
241  T factor = 1.0;
242  if (digits > 0)
243  {
244  for (int i = 0; i < digits; ++i)
245  factor *= 10.0;
246  }
247  else if (digits < 0)
248  {
249  for (int i = 0; i < -digits; ++i)
250  factor /= 10.0;
251  }
252 
253  return std::round(value * factor) / factor;
254  }
255 
273  template<typename T>
274  double percentOf(T value, T total, int digits)
275  {
276  if (value < 0) { throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Value must be non-negative", String(value)); }
277  if (total < 0) { throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Total must be non-negative", String(total)); }
278  if (total <= 0) // avoid float equality compare
279  {
280  return 0.0; // avoid division by zero
281  }
282  return roundTo(value * 100.0 / total, digits);
283  }
284 
290  inline bool approximatelyEqual(double a, double b, double tol)
291  {
292  return std::fabs(a - b) <= tol;
293  }
294 
303  template<typename T>
304  T gcd(T a, T b)
305  {
306  T c;
307  while (b != 0)
308  {
309  c = a % b;
310  a = b;
311  b = c;
312  }
313  return a;
314  }
315 
328  template<typename T>
329  T gcd(T a, T b, T& u1, T& u2)
330  {
331  u1 = 1;
332  u2 = 0;
333  T u3 = a;
334 
335  T v1 = 0;
336  T v2 = 1;
337  T v3 = b;
338 
339  while (v3 != 0)
340  {
341  T q = u3 / v3;
342  T t1 = u1 - v1 * q;
343  T t2 = u2 - v2 * q;
344  T t3 = u3 - v3 * q;
345 
346  u1 = v1;
347  u2 = v2;
348  u3 = v3;
349 
350  v1 = t1;
351  v2 = t2;
352  v3 = t3;
353  }
354 
355  return u3;
356  }
357 
367  template<typename T>
368  T getPPM(T mz_obs, T mz_ref)
369  {
370  return (mz_obs - mz_ref) / mz_ref * 1e6;
371  }
372 
382  template<typename T>
383  T getPPMAbs(T mz_obs, T mz_ref)
384  {
385  return std::fabs(getPPM(mz_obs, mz_ref));
386  }
387 
397  template<typename T>
398  T ppmToMass(T ppm, T mz_ref)
399  {
400  return (ppm / T(1e6)) * mz_ref;
401  }
402 
403  /*
404  @brief Compute the absolute mass diff in [Th], given a ppm value and a reference point.
405 
406  The returned mass diff is always positive!
407 
408  @param ppm Parts-per-million error
409  @param mz_ref Reference m/z
410  @return The absolute mass diff in [Th]
411  */
412  template<typename T>
413  T ppmToMassAbs(T ppm, T mz_ref)
414  {
415  return std::fabs(ppmToMass(ppm, mz_ref));
416  }
417 
431  inline std::pair<double, double> getTolWindow(double val, double tol, bool ppm)
432  {
433  double left, right;
434 
435  if (ppm)
436  {
437  left = val - val * tol * 1e-6;
438  right = val / (1.0 - tol * 1e-6);
439  }
440  else
441  {
442  left = val - tol;
443  right = val + tol;
444  }
445 
446  return std::make_pair(left, right);
447  }
448 
452  template<typename T1>
453  typename T1::value_type quantile(const T1& x, double q)
454  {
455  if (x.empty()) throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Quantile requested from empty container.");
456  if (q < 0.0) q = 0.;
457  if (q > 1.0) q = 1.;
458 
459  const auto n = x.size();
460  const auto id = std::max(0., n * q - 1); // -1 for c++ index starting at 0
461  const auto lo = floor(id);
462  const auto hi = ceil(id);
463  const auto qs = x[lo];
464  const auto h = (id - lo);
465 
466  return (1.0 - h) * qs + h * x[hi];
467  }
468 
469  // portable random shuffle
470  class OPENMS_DLLAPI RandomShuffler
471  {
472  public:
473  explicit RandomShuffler(int seed): rng_(boost::mt19937_64(seed))
474  {
475  }
476 
477  explicit RandomShuffler(const boost::mt19937_64& mt_rng): rng_(mt_rng)
478  {
479  }
480 
481  RandomShuffler() = default;
482  ~RandomShuffler() = default;
483 
484  boost::mt19937_64 rng_;
485  template<class RandomAccessIterator>
486  void portable_random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
487  {
488  for (auto i = (last - first) - 1; i > 0; --i) // OMS_CODING_TEST_EXCLUDE
489  {
490  boost::uniform_int<decltype(i)> d(0, i);
491  std::swap(first[i], first[d(rng_)]);
492  }
493  }
494 
495  void seed(uint64_t val)
496  {
497  rng_.seed(val);
498  }
499  };
500 } // namespace Math
501 } // namespace OpenMS
Exception indicating that an invalid parameter was handed over to an algorithm.
Definition: Exception.h:316
Invalid value exception.
Definition: Exception.h:305
Definition: MathFunctions.h:471
RandomShuffler(int seed)
Definition: MathFunctions.h:473
boost::mt19937_64 rng_
Definition: MathFunctions.h:484
void seed(uint64_t val)
Definition: MathFunctions.h:495
RandomShuffler(const boost::mt19937_64 &mt_rng)
Definition: MathFunctions.h:477
void portable_random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
Definition: MathFunctions.h:486
A more convenient string class.
Definition: String.h:34
unsigned int UInt
Unsigned integer type.
Definition: Types.h:64
#define OPENMS_PRECONDITION(condition, message)
Precondition macro.
Definition: openms/include/OpenMS/CONCEPT/Macros.h:94
T gcd(T a, T b)
Returns the greatest common divisor (gcd) of two numbers by applying the Euclidean algorithm.
Definition: MathFunctions.h:304
double intervalTransformation(double x, double left1, double right1, double left2, double right2)
transforms point x of interval [left1,right1] into interval [left2,right2]
Definition: MathFunctions.h:174
double linear2log(double x)
Transforms a number from linear to log10 scale. Avoids negative logarithms by adding 1.
Definition: MathFunctions.h:186
T round(T x)
Rounds the value.
Definition: MathFunctions.h:219
double log2linear(double x)
Transforms a number from log10 to to linear scale. Subtracts the 1 added by linear2log(double)
Definition: MathFunctions.h:198
double ceilDecimal(double x, int decPow)
rounds x up to the next decimal power 10 ^ decPow
Definition: MathFunctions.h:147
bool approximatelyEqual(double a, double b, double tol)
Returns if a is approximately equal b , allowing a tolerance of tol.
Definition: MathFunctions.h:290
double roundDecimal(double x, int decPow)
rounds x to the next decimal power 10 ^ decPow
Definition: MathFunctions.h:162
bool isOdd(UInt x)
Returns true if the given integer is odd.
Definition: MathFunctions.h:208
const double c
Definition: Constants.h:188
const double h
Definition: Constants.h:141
BinContainer createBins(double min, double max, uint32_t number_of_bins, double extend_margin=0)
Split a range [min,max] into number_of_bins (with optional overlap) and return the ranges of each bin...
Definition: MathFunctions.h:118
bool extendRange(T &min, T &max, const T &value)
Given an interval/range and a new value, extend the range to include the new value if needed.
Definition: MathFunctions.h:43
T getPPM(T mz_obs, T mz_ref)
Compute parts-per-million of two m/z values.
Definition: MathFunctions.h:368
T1::value_type quantile(const T1 &x, double q)
Returns the value of the q th quantile (0-1) in a sorted non-empty vector x.
Definition: MathFunctions.h:453
T roundTo(const T value, int digits)
Definition: MathFunctions.h:239
T ppmToMass(T ppm, T mz_ref)
Compute the mass diff in [Th], given a ppm value and a reference point.
Definition: MathFunctions.h:398
std::pair< double, double > zoomIn(const double left, const double right, const float factor, const float align)
Zoom into an interval [left, right], decreasing its width by factor (which must be in [0,...
Definition: MathFunctions.h:89
bool contains(T value, T min, T max)
Is a value contained in [min, max] ?
Definition: MathFunctions.h:64
T getPPMAbs(T mz_obs, T mz_ref)
Compute absolute parts-per-million of two m/z values.
Definition: MathFunctions.h:383
T ppmToMassAbs(T ppm, T mz_ref)
Definition: MathFunctions.h:413
double percentOf(T value, T total, int digits)
Definition: MathFunctions.h:274
std::pair< double, double > getTolWindow(double val, double tol, bool ppm)
Return tolerance window around val given tolerance tol.
Definition: MathFunctions.h:431
std::vector< RangeBase > BinContainer
Definition: MathFunctions.h:102
Main OpenMS namespace.
Definition: openswathalgo/include/OpenMS/OPENSWATHALGO/DATAACCESS/ISpectrumAccess.h:19
Base class for a simple range with minimum and maximum.
Definition: RangeManager.h:37