OpenMS  2.8.0
MathFunctions.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2021.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Timo Sachsenberg$
32 // $Authors: Marc Sturm $
33 // --------------------------------------------------------------------------
34 
35 #pragma once
36 
37 #include <OpenMS/CONCEPT/Types.h>
39 
40 #include <boost/random/mersenne_twister.hpp>
41 #include <boost/random/uniform_int.hpp>
42 #include <cmath>
43 #include <utility>
44 
45 namespace OpenMS
46 {
54  namespace Math
55  {
56 
65  template<typename T>
66  bool extendRange(T& min, T& max, const T& value)
67  {
68  if (value < min)
69  {
70  min = value;
71  return true;
72  }
73  if (value > max)
74  {
75  max = value;
76  return true;
77  }
78  return false;
79  }
80 
81 
93  inline static double ceilDecimal(double x, int decPow)
94  {
95  return (ceil(x / pow(10.0, decPow))) * pow(10.0, decPow); // decimal shift right, ceiling, decimal shift left
96  }
97 
108  inline static double roundDecimal(double x, int decPow)
109  {
110  if (x > 0)
111  return (floor(0.5 + x / pow(10.0, decPow))) * pow(10.0, decPow);
112 
113  return -((floor(0.5 + fabs(x) / pow(10.0, decPow))) * pow(10.0, decPow));
114  }
115 
121  inline static double intervalTransformation(double x, double left1, double right1, double left2, double right2)
122  {
123  return left2 + (x - left1) * (right2 - left2) / (right1 - left1);
124  }
125 
133  inline double linear2log(double x)
134  {
135  return log10(x + 1); //+1 to avoid negative logarithms
136  }
137 
145  inline double log2linear(double x)
146  {
147  return pow(10, x) - 1;
148  }
149 
155  inline bool isOdd(UInt x)
156  {
157  return (x & 1) != 0;
158  }
159 
165  template <typename T>
166  T round(T x)
167  {
168  if (x >= T(0))
169  {
170  return T(floor(x + T(0.5)));
171  }
172  else
173  {
174  return T(ceil(x - T(0.5)));
175  }
176  }
177 
183  inline static bool approximatelyEqual(double a, double b, double tol)
184  {
185  return std::fabs(a - b) <= tol;
186  }
187 
196  template <typename T>
197  T gcd(T a, T b)
198  {
199  T c;
200  while (b != 0)
201  {
202  c = a % b;
203  a = b;
204  b = c;
205  }
206  return a;
207  }
208 
221  template <typename T>
222  T gcd(T a, T b, T & u1, T & u2)
223  {
224  u1 = 1;
225  u2 = 0;
226  T u3 = a;
227 
228  T v1 = 0;
229  T v2 = 1;
230  T v3 = b;
231 
232  while (v3 != 0)
233  {
234  T q = u3 / v3;
235  T t1 = u1 - v1 * q;
236  T t2 = u2 - v2 * q;
237  T t3 = u3 - v3 * q;
238 
239  u1 = v1;
240  u2 = v2;
241  u3 = v3;
242 
243  v1 = t1;
244  v2 = t2;
245  v3 = t3;
246  }
247 
248  return u3;
249  }
250 
260  template <typename T>
261  T getPPM(T mz_obs, T mz_ref)
262  {
263  return (mz_obs - mz_ref) / mz_ref * 1e6;
264  }
265 
275  template <typename T>
276  T getPPMAbs(T mz_obs, T mz_ref)
277  {
278  return std::fabs(getPPM(mz_obs, mz_ref));
279  }
280 
290  template <typename T>
291  T ppmToMass(T ppm, T mz_ref)
292  {
293  return (ppm / 1e6) * mz_ref;
294  }
295 
296  /*
297  @brief Compute the absolute mass diff in [Th], given a ppm value and a reference point.
298 
299  The returned mass diff is always positive!
300 
301  @param ppm Parts-per-million error
302  @param mz_ref Reference m/z
303  @return The absolute mass diff in [Th]
304  */
305  template <typename T>
306  T ppmToMassAbs(T ppm, T mz_ref)
307  {
308  return std::fabs(ppmToMass(ppm, mz_ref));
309  }
310 
324  inline static std::pair<double, double> getTolWindow(double val, double tol, bool ppm)
325  {
326  double left, right;
327 
328  if (ppm)
329  {
330  left = val - val * tol * 1e-6;
331  right = val / (1.0 - tol * 1e-6);
332  }
333  else
334  {
335  left = val - tol;
336  right = val + tol;
337  }
338 
339  return std::make_pair(left, right);
340  }
341 
345  template <typename T1> typename T1::value_type quantile(const T1 &x, double q)
346  {
347  if (x.empty()) throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
348  "Quantile requested from empty container.");
349  if (q < 0.0) q = 0.;
350  if (q > 1.0) q = 1.;
351 
352  const auto n = x.size();
353  const auto id = std::max(0., n * q - 1); // -1 for c++ index starting at 0
354  const auto lo = floor(id);
355  const auto hi = ceil(id);
356  const auto qs = x[lo];
357  const auto h = (id - lo);
358 
359  return (1.0 - h) * qs + h * x[hi];
360  }
361 
362  // portable random shuffle
363  class OPENMS_DLLAPI RandomShuffler
364  {
365  public:
366  explicit RandomShuffler(int seed):
367  rng_(boost::mt19937_64(seed))
368  {}
369 
370  explicit RandomShuffler(const boost::mt19937_64& mt_rng):
371  rng_(mt_rng)
372  {}
373 
374  RandomShuffler() = default;
375  ~RandomShuffler() = default;
376 
377  boost::mt19937_64 rng_;
378  template <class RandomAccessIterator>
379  void portable_random_shuffle (RandomAccessIterator first, RandomAccessIterator last)
380  {
381  for (auto i = (last-first)-1; i > 0; --i) // OMS_CODING_TEST_EXCLUDE
382  {
383  boost::uniform_int<decltype(i)> d(0, i);
384  std::swap(first[i], first[d(rng_)]);
385  }
386  }
387 
388  void seed(uint64_t val)
389  {
390  rng_.seed(val);
391  }
392  };
393  } // namespace Math
394 } // namespace OpenMS
395 
Exception indicating that an invalid parameter was handed over to an algorithm.
Definition: Exception.h:341
Definition: MathFunctions.h:364
RandomShuffler(int seed)
Definition: MathFunctions.h:366
boost::mt19937_64 rng_
Definition: MathFunctions.h:377
void seed(uint64_t val)
Definition: MathFunctions.h:388
RandomShuffler(const boost::mt19937_64 &mt_rng)
Definition: MathFunctions.h:370
void portable_random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
Definition: MathFunctions.h:379
unsigned int UInt
Unsigned integer type.
Definition: Types.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:197
static double roundDecimal(double x, int decPow)
rounds x to the next decimal power 10 ^ decPow
Definition: MathFunctions.h:108
double linear2log(double x)
Transforms a number from linear to log10 scale. Avoids negative logarithms by adding 1.
Definition: MathFunctions.h:133
T round(T x)
Rounds the value.
Definition: MathFunctions.h:166
double log2linear(double x)
Transforms a number from log10 to to linear scale. Subtracts the 1 added by linear2log(double)
Definition: MathFunctions.h:145
static 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:121
static double ceilDecimal(double x, int decPow)
rounds x up to the next decimal power 10 ^ decPow
Definition: MathFunctions.h:93
static bool approximatelyEqual(double a, double b, double tol)
Returns if a is approximately equal b , allowing a tolerance of tol.
Definition: MathFunctions.h:183
bool isOdd(UInt x)
Returns true if the given integer is odd.
Definition: MathFunctions.h:155
const double c
Definition: Constants.h:209
const double h
Definition: Constants.h:162
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:66
T getPPM(T mz_obs, T mz_ref)
Compute parts-per-million of two m/z values.
Definition: MathFunctions.h:261
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:345
static std::pair< double, double > getTolWindow(double val, double tol, bool ppm)
Return tolerance window around val given tolerance tol.
Definition: MathFunctions.h:324
T ppmToMass(T ppm, T mz_ref)
Compute the mass diff in [Th], given a ppm value and a reference point.
Definition: MathFunctions.h:291
T getPPMAbs(T mz_obs, T mz_ref)
Compute absolute parts-per-million of two m/z values.
Definition: MathFunctions.h:276
T ppmToMassAbs(T ppm, T mz_ref)
Definition: MathFunctions.h:306
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47