Sometimes a value gets negative zero causing a NaN on the sqrt.

Closes #61

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;
   }