about summary refs log tree commit diff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/debug.cpp6
-rw-r--r--src/io.cpp64
-rw-r--r--src/lmm.cpp6
-rw-r--r--src/mathfunc.cpp2
-rw-r--r--src/mathfunc.h3
-rwxr-xr-xtest/dev_test_suite.sh3
6 files changed, 46 insertions, 38 deletions
diff --git a/src/debug.cpp b/src/debug.cpp
index a00190d..cb98f27 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -20,7 +20,7 @@
 
 static bool debug_mode     = false;
 static bool debug_check    = true;  // check data/algorithms
-static bool debug_strict   = false; // fail on error
+static bool debug_strict   = false; // fail on error, more rigorous checks
 static bool debug_quiet    = false;
 static uint debug_issue    = 0;     // track github issues
 static bool debug_legacy   = false; // legacy mode
@@ -64,7 +64,7 @@ int gsl_matrix_safe_memcpy (gsl_matrix *dest, const gsl_matrix *src) {
 
 void do_gsl_matrix_safe_free (gsl_matrix *m, const char *__pretty_function, const char *__file, int __line) {
   enforce(m);
-  if (is_check_mode() && is_debug_mode()) {
+  if (is_strict_mode() && is_check_mode() && is_debug_mode()) {
     bool has_NaN = has_nan(m);
     bool has_Inf = has_inf(m);
     if (has_NaN || has_Inf) {
@@ -89,7 +89,7 @@ int gsl_vector_safe_memcpy (gsl_vector *dest, const gsl_vector *src) {
 
 void do_gsl_vector_safe_free (gsl_vector *v, const char *__pretty_function, const char *__file, int __line) {
   enforce(v);
-  if (is_check_mode() && is_debug_mode()) {
+  if (is_strict_mode() && is_check_mode() && is_debug_mode()) {
     bool has_NaN = has_nan(v);
     bool has_Inf = has_inf(v);
     if (has_NaN || has_Inf) {
diff --git a/src/io.cpp b/src/io.cpp
index ef6c26a..fb0f0e7 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -3265,8 +3265,7 @@ void ReadFile_beta(const string &file_beta,
   string type;
 
   string rs, chr, a1, a0, pos, cm;
-  double z = 0, beta = 0, se_beta = 0, chisq = 0, pvalue = 0, zsquare = 0,
-         af = 0, var_x = 0;
+  double z = 0, beta = 0, se_beta = 0, pvalue = 0, zsquare = 0; // af = 0;
   size_t n_total = 0, n_mis = 0, n_obs = 0, n_case = 0, n_control = 0;
 
   // Read header.
@@ -3298,15 +3297,15 @@ void ReadFile_beta(const string &file_beta,
     z = 0;
     beta = 0;
     se_beta = 0;
-    chisq = 0;
+    auto chisq = 0.0;
     pvalue = 0;
     n_total = 0;
     n_mis = 0;
     n_obs = 0;
     n_case = 0;
     n_control = 0;
-    af = 0;
-    var_x = 0;
+    // af = 0;
+    // auto var_x = 0.0;
     for (size_t i = 0; i < header.coln; i++) {
       enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
@@ -3338,7 +3337,7 @@ void ReadFile_beta(const string &file_beta,
         se_beta = atof(ch_ptr);
       }
       if (header.chisq_col != 0 && header.chisq_col == i + 1) {
-        chisq = atof(ch_ptr);
+         chisq = atof(ch_ptr);
       }
       if (header.p_col != 0 && header.p_col == i + 1) {
         pvalue = atof(ch_ptr);
@@ -3359,12 +3358,12 @@ void ReadFile_beta(const string &file_beta,
       if (header.ncontrol_col != 0 && header.ncontrol_col == i + 1) {
         n_control = atoi(ch_ptr);
       }
-      if (header.af_col != 0 && header.af_col == i + 1) {
-        af = atof(ch_ptr);
-      }
-      if (header.var_col != 0 && header.var_col == i + 1) {
-        var_x = atof(ch_ptr);
-      }
+      // if (header.af_col != 0 && header.af_col == i + 1) {
+      //   af = atof(ch_ptr);
+      // }
+      // if (header.var_col != 0 && header.var_col == i + 1) {
+      //   var_x = atof(ch_ptr);
+      // }
 
       ch_ptr = strtok(NULL, " , \t");
     }
@@ -3397,9 +3396,9 @@ void ReadFile_beta(const string &file_beta,
     }
 
     // Obtain var_x.
-    if (header.var_col == 0 && header.af_col != 0) {
-      var_x = 2.0 * af * (1.0 - af);
-    }
+    // if (header.var_col == 0 && header.af_col != 0) {
+    //   var_x = 2.0 * af * (1.0 - af);
+    // }
 
     // If the SNP is also present in cor file, then do calculations.
     if ((mapRS2wA.size() == 0 || mapRS2wA.count(rs) != 0) &&
@@ -3448,7 +3447,7 @@ void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
   string type;
 
   string rs, chr, a1, a0, pos, cm;
-  double z = 0, beta = 0, se_beta = 0, chisq = 0, pvalue = 0, af = 0, var_x = 0;
+  double z = 0, beta = 0, se_beta = 0; // pvalue = 0, chisq=0, af = 0 , var_x = 0;
   size_t n_total = 0, n_mis = 0, n_obs = 0, n_case = 0, n_control = 0;
   size_t ni_total = 0, ns_total = 0, ns_test = 0;
 
@@ -3480,15 +3479,15 @@ void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
     z = 0;
     beta = 0;
     se_beta = 0;
-    chisq = 0;
-    pvalue = 0;
+    // chisq = 0;
+    // pvalue = 0;
     n_total = 0;
     n_mis = 0;
     n_obs = 0;
     n_case = 0;
     n_control = 0;
-    af = 0;
-    var_x = 0;
+    // af = 0;
+    // double var_x = 0;
     for (size_t i = 0; i < header.coln; i++) {
       enforce(ch_ptr);
       if (header.rs_col != 0 && header.rs_col == i + 1) {
@@ -3519,12 +3518,12 @@ void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
       if (header.sebeta_col != 0 && header.sebeta_col == i + 1) {
         se_beta = atof(ch_ptr);
       }
-      if (header.chisq_col != 0 && header.chisq_col == i + 1) {
-        chisq = atof(ch_ptr);
-      }
-      if (header.p_col != 0 && header.p_col == i + 1) {
-        pvalue = atof(ch_ptr);
-      }
+      // if (header.chisq_col != 0 && header.chisq_col == i + 1) {
+      //   chisq = atof(ch_ptr);
+      // }
+      // if (header.p_col != 0 && header.p_col == i + 1) {
+      //   pvalue = atof(ch_ptr);
+      // }
 
       if (header.n_col != 0 && header.n_col == i + 1) {
         n_total = atoi(ch_ptr);
@@ -3542,12 +3541,13 @@ void ReadFile_beta(const string &file_beta, const map<string, double> &mapRS2wA,
         n_control = atoi(ch_ptr);
       }
 
-      if (header.af_col != 0 && header.af_col == i + 1) {
-        af = atof(ch_ptr);
-      }
-      if (header.var_col != 0 && header.var_col == i + 1) {
-        var_x = atof(ch_ptr);
-      }
+      // if (header.af_col != 0 && header.af_col == i + 1) {
+      //   af = atof(ch_ptr);
+      // }
+
+      // if (header.var_col != 0 && header.var_col == i + 1) {
+      //   var_x = atof(ch_ptr);
+      // }
 
       ch_ptr = strtok(NULL, " , \t");
     }
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 1a2614b..99a8436 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -45,6 +45,7 @@
 #include "fastblas.h"
 #include "lapack.h"
 #include "lmm.h"
+#include "mathfunc.h"
 
 using namespace std;
 
@@ -228,6 +229,7 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
               gsl_matrix_const_column(Uab, index_ab);
           gsl_blas_ddot(Hi_eval, &Uab_col.vector, &p_ab);
           if (e_mode != 0) {
+            assert(false);
             p_ab = gsl_vector_get(ab, index_ab) - p_ab;
           }
           gsl_matrix_set(Pab, 0, index_ab, p_ab);
@@ -389,8 +391,10 @@ double LogL_f(double l, void *params) {
 
   index_yy = GetabIndex(n_cvt + 2, n_cvt + 2, n_cvt);
   double P_yy = gsl_matrix_get(Pab, nc_total, index_yy);
-  f = c - 0.5 * logdet_h - 0.5 * (double)ni_test * log(P_yy);
 
+  assert(!is_nan(P_yy));
+  f = c - 0.5 * logdet_h - 0.5 * (double)ni_test * log(P_yy);
+  assert(!is_nan(f));
   gsl_matrix_safe_free(Pab); // FIXME
   gsl_vector_safe_free(Hi_eval);
   gsl_vector_safe_free(v_temp);
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index 908a52d..d3804dd 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -58,7 +58,7 @@ using namespace std;
 
 bool has_nan(const vector<double> v) {
   for (const auto& e: v) {
-    if (std::isnan(e))
+    if (is_nan(e))
       return true;
   }
   return false;
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 42e76e5..0d791a1 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -28,6 +28,9 @@
 
 using namespace std;
 
+inline bool is_nan(double f) {
+  return (std::isnan(f));
+}
 
 bool has_nan(const vector<double> v);
 bool has_nan(const gsl_vector *v);
diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh
index 851d67f..0d3d8a0 100755
--- a/test/dev_test_suite.sh
+++ b/test/dev_test_suite.sh
@@ -1,7 +1,8 @@
 #!/usr/bin/env bash
 
 gemma=../bin/gemma
-gemmaopts=-debug
+# gemmaopts="-debug -strict"
+gemmaopts="-debug"
 
 # Related to https://github.com/genetics-statistics/GEMMA/issues/78
 testBXDStandardRelatednessMatrixKSingularError() {
generated by cgit-pink 1.4.1 (git 2.36.1) at 2025-10-07 21:44:05 +0000