IsoSpec 2.2.1
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
34namespace IsoSpec
35{
36
37extern double* g_lfact_table;
38
39static 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
53const double pi = 3.14159265358979323846264338328;
54const double logpi = 1.144729885849400174143427351353058711647294812915311571513623071472137769884826079783623270275489708;
55
56double NormalCDFInverse(double p);
57double NormalCDFInverse(double p, double mean, double stdev);
58double NormalCDF(double x, double mean, double stdev);
59double 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.
62double LowerIncompleteGamma2(int a, double x);
63
64// Returns y such that LowerIncompleteGamma2(a, y) == x. Approximately.
65double InverseLowerIncompleteGamma2(int a, double x);
66
67// Computes the inverse Cumulative Distribution Funcion of the Chi-Square distribution with k degrees of freedom
68inline double InverseChiSquareCDF2(int k, double x)
69{
70 return InverseLowerIncompleteGamma2(k, x*tgamma(static_cast<double>(k)/2.0)) * 2.0;
71}
72
73extern std::mt19937 random_gen;
74extern std::uniform_real_distribution<double> stdunif;
75
76inline 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
82size_t rdvariate_binom(size_t tries, double succ_prob, std::mt19937& rgen = random_gen);
83
84
85
86
87} // namespace IsoSpec