IsoSpec  2.1.2
isoMath.h
1 /*
2  * Copyright (C) 2015-2020 Mateusz Łącki and Michał Startek.
3  *
4  * This file is part of IsoSpec.
5  *
6  * IsoSpec is free software: you can redistribute it and/or modify
7  * it under the terms of the Simplified ("2-clause") BSD licence.
8  *
9  * IsoSpec is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
12  *
13  * You should have received a copy of the Simplified BSD Licence
14  * along with IsoSpec. If not, see <https://opensource.org/licenses/BSD-2-Clause>.
15  */
16 
17 #pragma once
18 
19 #include <cmath>
20 #include <random>
21 
22 #if !defined(ISOSPEC_G_FACT_TABLE_SIZE)
23 // 10M should be enough for anyone, right?
24 // Actually, yes. If anyone tries to input a molecule that has more than 10M atoms,
25 // he deserves to get an exception thrown in his face. OpenMS guys don't want to alloc
26 // a table of 10M to memoize the necessary values though, use something smaller for them.
27  #if ISOSPEC_BUILDING_OPENMS
28  #define ISOSPEC_G_FACT_TABLE_SIZE 1024
29  #else
30  #define ISOSPEC_G_FACT_TABLE_SIZE 1024*1024*10
31  #endif
32 #endif
33 
34 namespace IsoSpec
35 {
36 
37 extern double* g_lfact_table;
38 
39 static inline double minuslogFactorial(int n)
40 {
41  if (n < 2)
42  return 0.0;
43  #if ISOSPEC_BUILDING_OPENMS
44  if (n >= ISOSPEC_G_FACT_TABLE_SIZE)
45  return -lgamma(n+1);
46  #endif
47  if (g_lfact_table[n] == 0.0)
48  g_lfact_table[n] = -lgamma(n+1);
49 
50  return g_lfact_table[n];
51 }
52 
53 const double pi = 3.14159265358979323846264338328;
54 const double logpi = 1.144729885849400174143427351353058711647294812915311571513623071472137769884826079783623270275489708;
55 
56 double NormalCDFInverse(double p);
57 double NormalCDFInverse(double p, double mean, double stdev);
58 double NormalCDF(double x, double mean, double stdev);
59 double NormalPDF(double x, double mean = 0.0, double stdev = 1.0);
60 
61 // Returns lower incomplete gamma function of a/2, x, where a is int and > 0.
62 double LowerIncompleteGamma2(int a, double x);
63 
64 // Returns y such that LowerIncompleteGamma2(a, y) == x. Approximately.
65 double InverseLowerIncompleteGamma2(int a, double x);
66 
67 // Computes the inverse Cumulative Distribution Funcion of the Chi-Square distribution with k degrees of freedom
68 inline double InverseChiSquareCDF2(int k, double x)
69 {
70  return InverseLowerIncompleteGamma2(k, x*tgamma(static_cast<double>(k)/2.0)) * 2.0;
71 }
72 
73 extern std::mt19937 random_gen;
74 extern std::uniform_real_distribution<double> stdunif;
75 
76 inline double rdvariate_beta_1_b(double b, std::mt19937& rgen = random_gen)
77 {
78  return 1.0 - pow(stdunif(rgen), 1.0/b);
79 }
80 
81 
82 size_t rdvariate_binom(size_t tries, double succ_prob, std::mt19937& rgen = random_gen);
83 
84 
85 
86 
87 } // namespace IsoSpec
IsoSpec
Definition: allocator.cpp:20