1 /*
2 Genome-wide Efficient Mixed Model Association (GEMMA)
3 Copyright © 2011-2017, Xiang Zhou
4 Copyright © 2017, Peter Carbonetto
5 Copyright © 2017, Pjotr Prins
6
7 This program is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #ifndef __MATHFUNC_H__
22 #define __MATHFUNC_H__
23
24 // #include "Eigen/Dense"
25 #include "gsl/gsl_matrix.h"
26 #include "gsl/gsl_vector.h"
27
28 #define CONDITIONED_MAXRATIO 2e+6 // based on http://mathworld.wolfram.com/ConditionNumber.html
29 #define EIGEN_MINVALUE 1e-10
30
31 using namespace std;
32
is_nan(double f)33 inline bool is_nan(double f) {
34 return (isnan(f));
35 }
36
37 #define is_inf(d) gsl_isinf(d)
38 #define is_nan(d) gsl_isnan(d)
39
40 bool has_nan(const vector<double> v);
41 bool has_nan(const gsl_vector *v);
42 bool has_inf(const gsl_vector *v);
43 bool has_nan(const gsl_matrix *m);
44 bool has_inf(const gsl_matrix *m);
45
46 bool is_integer(const std::string & s);
47
48 double safe_log(const double d);
49 double safe_sqrt(const double d);
50
51 double VectorVar(const gsl_vector *v);
52 void CenterMatrix(gsl_matrix *G);
53 void CenterMatrix(gsl_matrix *G, const gsl_vector *w);
54 void CenterMatrix(gsl_matrix *G, const gsl_matrix *W);
55 void StandardizeMatrix(gsl_matrix *G);
56 double ScaleMatrix(gsl_matrix *G);
57 bool has_negative_values_but_one(const gsl_vector *v);
58 uint count_abs_small_values(const gsl_vector *v, double min);
59 bool isMatrixPositiveDefinite(const gsl_matrix *G);
60 bool isMatrixIllConditioned(const gsl_vector *eigenvalues, double max_ratio=CONDITIONED_MAXRATIO);
61 bool isMatrixSymmetric(const gsl_matrix *G);
62 gsl_vector *getEigenValues(const gsl_matrix *G);
63 double sum(const double *m, size_t rows, size_t cols);
64 double SumVector(const gsl_vector *v);
65 double CenterVector(gsl_vector *y);
66 void CenterVector(gsl_vector *y, const gsl_matrix *W);
67 void StandardizeVector(gsl_vector *y);
68 void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX);
69 void CalcUtX(const gsl_matrix *U, const gsl_matrix *X, gsl_matrix *UtX);
70 void CalcUtX(const gsl_matrix *U, const gsl_vector *x, gsl_vector *Utx);
71 double CalcHWE(const size_t n_hom1, const size_t n_hom2, const size_t n_ab);
72 void Kronecker(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H);
73 void KroneckerSym(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H);
74
75 double UcharToDouble02(const unsigned char c);
76 unsigned char Double02ToUchar(const double dosage);
77 // void uchar_matrix_get_row(const vector<vector<unsigned char>> &X,
78 // const size_t i_row, Eigen::VectorXd &x_row);
79
80 #endif
81