aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/gemma.cpp12
-rw-r--r--src/lmm.cpp10
-rw-r--r--src/mathfunc.cpp13
-rw-r--r--src/mathfunc.h1
-rw-r--r--src/mvlmm.cpp47
5 files changed, 51 insertions, 32 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp
index cb01ee0..f928977 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2595,7 +2595,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
if (wi <= 0 || wj <= 0) {
d = 0;
} else {
- d /= sqrt(wi * wj);
+ d /= safe_sqrt(wi * wj);
}
gsl_matrix_set(G, i, j, d);
if (j != i) {
@@ -2623,7 +2623,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
if (wi <= 0) {
wi = 0;
} else {
- wi = sqrt(wi);
+ wi = safe_sqrt(wi);
}
gsl_vector_view Urow = gsl_matrix_row(U, i);
gsl_vector_scale(&Urow.vector, wi);
@@ -3418,7 +3418,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
for (size_t j = 0; j <= i; j++) {
c = (2 * cPar.n_ph - min(i, j) + 1) * min(i, j) / 2 + max(i, j) -
min(i, j);
- outfile << tab(j) << sqrt(cPar.VVg_remle_null[c]);
+ outfile << tab(j) << safe_sqrt(cPar.VVg_remle_null[c]);
}
outfile << endl;
}
@@ -3436,7 +3436,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
for (size_t j = 0; j <= i; j++) {
c = (2 * cPar.n_ph - min(i, j) + 1) * min(i, j) / 2 + max(i, j) -
min(i, j);
- outfile << tab(j) << sqrt(cPar.VVe_remle_null[c]);
+ outfile << tab(j) << safe_sqrt(cPar.VVe_remle_null[c]);
}
outfile << endl;
}
@@ -3455,7 +3455,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
for (size_t j = 0; j <= i; j++) {
c = (2 * cPar.n_ph - min(i, j) + 1) * min(i, j) / 2 + max(i, j) -
min(i, j);
- outfile << tab(j) << sqrt(cPar.VVg_mle_null[c]);
+ outfile << tab(j) << safe_sqrt(cPar.VVg_mle_null[c]);
}
outfile << endl;
}
@@ -3473,7 +3473,7 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) {
for (size_t j = 0; j <= i; j++) {
c = (2 * cPar.n_ph - min(i, j) + 1) * min(i, j) / 2 + max(i, j) -
min(i, j);
- outfile << tab(j) << sqrt(cPar.VVe_mle_null[c]);
+ outfile << tab(j) << safe_sqrt(cPar.VVe_mle_null[c]);
}
outfile << endl;
}
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 80372ee..dbf1ad5 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1144,7 +1144,7 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM &params, double &beta,
beta = P_xy / P_xx;
double tau = (double)df / Px_yy;
- se = sqrt(1.0 / (tau * P_xx));
+ se = safe_sqrt(1.0 / (tau * P_xx));
p_wald = gsl_cdf_fdist_Q((P_yy - Px_yy) * tau, 1.0, df);
gsl_matrix_free(Pab);
@@ -1186,7 +1186,7 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM &params, double &beta,
beta = P_xy / P_xx;
double tau = (double)df / Px_yy;
- se = sqrt(1.0 / (tau * P_xx));
+ se = safe_sqrt(1.0 / (tau * P_xx));
p_score =
gsl_cdf_fdist_Q((double)ni_test * P_xy * P_xy / (P_yy * P_xx), 1.0, df);
@@ -1201,7 +1201,7 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
size_t index_ab;
size_t n_cvt = UtW->size2;
- debug_msg("entering");
+ // debug_msg("entering");
gsl_vector *u_a = gsl_vector_safe_alloc(Uty->size);
@@ -2018,7 +2018,7 @@ void CalcPve(const gsl_vector *eval, const gsl_matrix *UtW,
FUNC_PARAM param0 = {true, ni_test, n_cvt, eval, Uab, ab, 0};
- double se = sqrt(-1.0 / LogRL_dev2(lambda, &param0));
+ double se = safe_sqrt(-1.0 / LogRL_dev2(lambda, &param0));
pve = trace_G * lambda / (trace_G * lambda + 1.0);
pve_se = trace_G / ((trace_G * lambda + 1.0) * (trace_G * lambda + 1.0)) * se;
@@ -2085,7 +2085,7 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
// Obtain se_beta.
for (size_t i = 0; i < Vbeta->size1; i++) {
- gsl_vector_set(se_beta, i, sqrt(gsl_matrix_get(Vbeta, i, i)));
+ gsl_vector_set(se_beta, i, safe_sqrt(gsl_matrix_get(Vbeta, i, i)));
}
gsl_matrix_safe_free(Uab);
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index fdd6a58..542093e 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -104,11 +104,22 @@ bool is_float(const std::string & s){
}
double safe_log(const double d) {
- if (!is_legacy_mode())
+ if (!is_legacy_mode() && !is_check_mode())
enforce_msg(d > 0.0, (std::string("Trying to take the log of ") + std::to_string(d)).c_str());
return log(d);
}
+double safe_sqrt(const double d) {
+ double d1 = d;
+ if (fabs(d < 0.001))
+ d1 = fabs(d);
+ if (!is_legacy_mode() && !is_check_mode())
+ enforce_msg(d1 >= 0.0, (std::string("Trying to take the sqrt of ") + std::to_string(d)).c_str());
+ if (d1 < 0.0 )
+ return nan("");
+ return sqrt(d1);
+}
+
// calculate variance of a vector
double VectorVar(const gsl_vector *v) {
double d, m = 0.0, m2 = 0.0;
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 5b1e2ec..e786e87 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -46,6 +46,7 @@ bool has_inf(const gsl_matrix *m);
bool is_integer(const std::string & s);
double safe_log(const double d);
+double safe_sqrt(const double d);
double VectorVar(const gsl_vector *v);
void CenterMatrix(gsl_matrix *G);
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index faea21e..a3c7c85 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -237,7 +237,7 @@ double EigenProc(const gsl_matrix *V_g, const gsl_matrix *V_e, gsl_vector *D_l,
logdet_Ve += safe_log(d);
gsl_vector_view U_col = gsl_matrix_column(U_l, i);
- d = sqrt(d);
+ d = safe_sqrt(d);
gsl_blas_dsyr(CblasUpper, d, &U_col.vector, V_e_h);
d = 1.0 / d;
gsl_blas_dsyr(CblasUpper, d, &U_col.vector, V_e_hi);
@@ -326,7 +326,7 @@ double CalcQi(const gsl_vector *eval, const gsl_vector *D_l,
void CalcXHiY(const gsl_vector *eval, const gsl_vector *D_l,
const gsl_matrix *X, const gsl_matrix *UltVehiY,
gsl_vector *xHiy) {
- debug_msg("enter");
+ // debug_msg("enter");
size_t n_size = eval->size, c_size = X->size1, d_size = D_l->size;
// gsl_vector_set_zero(xHiy);
@@ -907,7 +907,7 @@ void MphCalcBeta(const gsl_vector *eval, const gsl_matrix *W,
for (size_t j = 0; j < B->size2; j++) {
for (size_t i = 0; i < B->size1; i++) {
gsl_matrix_set(B, i, j, gsl_vector_get(beta, j * d_size + i));
- gsl_matrix_set(se_B, i, j, sqrt(gsl_matrix_get(Vbeta, j * d_size + i,
+ gsl_matrix_set(se_B, i, j, safe_sqrt(gsl_matrix_get(Vbeta, j * d_size + i,
j * d_size + i)));
}
}
@@ -2948,7 +2948,7 @@ double PCRT(const size_t mode, const size_t d_size, const double p_value,
if (mode == 1) {
double a = crt_c / (2.0 * q * (q + 2.0));
double b = 1.0 + (crt_a + crt_b) / (2.0 * q);
- chisq_crt = (-1.0 * b + sqrt(b * b + 4.0 * a * chisq)) / (2.0 * a);
+ chisq_crt = (-1.0 * b + safe_sqrt(b * b + 4.0 * a * chisq)) / (2.0 * a);
} else if (mode == 2) {
chisq_crt = chisq / (1.0 + crt_a / (2.0 * q));
} else {
@@ -3090,7 +3090,7 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -3105,7 +3105,7 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
}
cout << endl;
}
@@ -3153,7 +3153,7 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -3168,7 +3168,7 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
}
cout << endl;
}
@@ -3529,7 +3529,7 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -3544,7 +3544,11 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ auto v = gsl_matrix_get(Hessian, c + v_size, c + v_size);
+ if (is_strict_mode())
+ enforce_msg(v >= 0,"se(Ve) is not valid");
+
+ cout << safe_sqrt(v) << "\t";
}
cout << endl;
}
@@ -3592,7 +3596,7 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -3607,7 +3611,10 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ auto v = gsl_matrix_get(Hessian, c + v_size, c + v_size);
+ if (is_strict_mode())
+ enforce_msg(v >= 0,"se(Ve) is not valid");
+ cout << safe_sqrt(v) << "\t";
}
cout << endl;
}
@@ -4076,7 +4083,7 @@ void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -4091,7 +4098,7 @@ void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
}
cout << endl;
}
@@ -4139,7 +4146,7 @@ void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -4154,7 +4161,7 @@ void MVLMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
}
cout << endl;
}
@@ -4521,7 +4528,7 @@ void MVLMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -4536,7 +4543,7 @@ void MVLMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
}
cout << endl;
}
@@ -4584,7 +4591,7 @@ void MVLMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c, c)) << "\t";
}
cout << endl;
}
@@ -4599,7 +4606,7 @@ void MVLMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval,
for (size_t i = 0; i < d_size; i++) {
for (size_t j = 0; j <= i; j++) {
c = GetIndex(i, j, d_size);
- cout << sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
+ cout << safe_sqrt(gsl_matrix_get(Hessian, c + v_size, c + v_size)) << "\t";
}
cout << endl;
}
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-18 19:27:29 +0000