about summary refs log tree commit diff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/debug.cpp19
-rw-r--r--src/debug.h2
-rw-r--r--src/gemma.cpp4
-rw-r--r--src/lmm.cpp48
-rw-r--r--src/param.cpp8
5 files changed, 52 insertions, 29 deletions
diff --git a/src/debug.cpp b/src/debug.cpp
index 3efcce6..bc40b44 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -38,15 +38,17 @@
 #include "debug.h"
 #include "mathfunc.h"
 
-static bool debug_mode     = false;
-static bool debug_check    = true;  // check data/algorithms
-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
+static bool debug_mode      = false;
+static bool debug_check     = true;  // check data/algorithms
+static bool debug_fpe_check = true;  // check floating point errors (intel hardware)
+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
 
 void debug_set_debug_mode(bool setting) { debug_mode = setting; }
 void debug_set_no_check_mode(bool setting) {debug_check = !setting; }
+void debug_set_no_fpe_check_mode(bool setting) {debug_fpe_check = !setting; }
 void debug_set_strict_mode(bool setting) { debug_strict = setting; }
 void debug_set_quiet_mode(bool setting) { debug_quiet = setting; }
 void debug_set_issue(uint issue) { debug_issue = issue; }
@@ -55,6 +57,7 @@ void debug_set_legacy_mode(bool setting) { debug_legacy = setting; }
 bool is_debug_mode() { return debug_mode; };
 bool is_no_check_mode() { return !debug_check; };
 bool is_check_mode() { return debug_check; };
+bool is_fpe_check_mode() { return debug_fpe_check; };
 bool is_strict_mode() { return debug_strict; };
 bool is_quiet_mode() { return debug_quiet; };
 bool is_issue(uint issue) { return issue == debug_issue; };
@@ -128,7 +131,7 @@ inline int fedisableexcept(unsigned int excepts)
 #endif
 
 void enable_segfpe() {
-  if (is_legacy_mode()) return;
+  if (!is_fpe_check_mode() || is_legacy_mode()) return;
   #ifdef __GNUC__
     #if defined(__x86_64__)
        debug_msg("enable segfpe hardware floating point error detection");
@@ -139,7 +142,7 @@ void enable_segfpe() {
 }
 
 void disable_segfpe() {
-  if (is_legacy_mode()) return;
+  if (!is_fpe_check_mode() || is_legacy_mode()) return;
   #ifdef __GNUC__
     #if defined(__x86_64__)
       debug_msg("disable segfpe");
diff --git a/src/debug.h b/src/debug.h
index a4be9f5..3133963 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -32,6 +32,7 @@ void gemma_gsl_error_handler (const char * reason,
 
 void debug_set_debug_mode(bool setting);
 void debug_set_no_check_mode(bool setting);
+void debug_set_no_fpe_check_mode(bool setting);
 void debug_set_strict_mode(bool setting);
 void debug_set_quiet_mode(bool setting);
 void debug_set_issue(uint issue);
@@ -40,6 +41,7 @@ void debug_set_legacy_mode(bool setting);
 bool is_debug_mode();
 bool is_no_check_mode();
 bool is_check_mode();
+bool is_fpe_check_mode();
 bool is_strict_mode();
 bool is_quiet_mode();
 bool is_issue(uint issue);
diff --git a/src/gemma.cpp b/src/gemma.cpp
index e1d76d4..7d12055 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -733,6 +733,7 @@ void GEMMA::PrintHelp(size_t option) {
   if (option == 14) {
     cout << " DEBUG OPTIONS" << endl;
     cout << " -no-check                disable checks" << endl;
+    cout << " -no-fpe-check            disable hardware floating point checking" << endl;
     cout << " -strict                  strict mode will stop when there is a problem" << endl;
     cout << " -silence                 silent terminal display" << endl;
     cout << " -debug                   debug output" << endl;
@@ -1616,6 +1617,9 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
     } else if (strcmp(argv[i], "-no-check") == 0) {
       // cPar.mode_check = false;
       debug_set_no_check_mode(true);
+    } else if (strcmp(argv[i], "-no-fpe-check") == 0) {
+      // cPar.mode_check = false;
+      debug_set_no_fpe_check_mode(true);
     } else if (strcmp(argv[i], "-strict") == 0) {
       // cPar.mode_strict = true;
       debug_set_strict_mode(true);
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 3f222d1..0fa8ea5 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -330,10 +330,11 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
           double ps_ww = gsl_matrix_safe_get(Pab, p - 1, index_ww);
 
           cout << "unsafe " << p-1 << "," << index_ww << ":" << gsl_matrix_get(Pab,p-1,index_ww) << endl;
-          assert(ps_ww != 0.0);
-          p_ab = ps_ab - ps_aw * ps_bw / ps_ww;
+          // if (is_check_mode() || is_debug_mode()) assert(ps_ww != 0.0);
+          if (ps_ww != 0)
+            p_ab = ps_ab - ps_aw * ps_bw / ps_ww;
 
-          cout << p << "r," << index_ab << ":" << p_ab << endl;
+          cout << "set " << p << "r," << index_ab << "c: " << p_ab << endl;
           gsl_matrix_set(Pab, p, index_ab, p_ab);
         }
       }
@@ -384,10 +385,13 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode,
           ps2_bw = gsl_matrix_safe_get(PPab, p - 1, index_bw);
           ps2_ww = gsl_matrix_safe_get(PPab, p - 1, index_ww);
 
-          assert(ps_ww != 0.0);
-          p2_ab = ps2_ab + ps_aw * ps_bw * ps2_ww / (ps_ww * ps_ww);
-          p2_ab -= (ps_aw * ps2_bw + ps_bw * ps2_aw) / ps_ww;
+          // if (is_check_mode() || is_debug_mode()) assert(ps_ww != 0.0);
+          if (ps_ww != 0) {
+            p2_ab = ps2_ab + ps_aw * ps_bw * ps2_ww / (ps_ww * ps_ww);
+            p2_ab -= (ps_aw * ps2_bw + ps_bw * ps2_aw) / ps_ww;
+          }
           gsl_matrix_set(PPab, p, index_ab, p2_ab);
+
         }
       }
     }
@@ -441,14 +445,15 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode,
           ps3_bw = gsl_matrix_safe_get(PPPab, p - 1, index_bw);
           ps3_ww = gsl_matrix_safe_get(PPPab, p - 1, index_ww);
 
-          assert(ps_ww != 0.0);
-          p3_ab = ps3_ab -
-                  ps_aw * ps_bw * ps2_ww * ps2_ww / (ps_ww * ps_ww * ps_ww);
-          p3_ab -= (ps_aw * ps3_bw + ps_bw * ps3_aw + ps2_aw * ps2_bw) / ps_ww;
-          p3_ab += (ps_aw * ps2_bw * ps2_ww + ps_bw * ps2_aw * ps2_ww +
-                    ps_aw * ps_bw * ps3_ww) /
-                   (ps_ww * ps_ww);
-
+          // if (is_check_mode() || is_debug_mode()) assert(ps_ww != 0.0);
+          if (ps_ww != 0) {
+            p3_ab = ps3_ab -
+              ps_aw * ps_bw * ps2_ww * ps2_ww / (ps_ww * ps_ww * ps_ww);
+            p3_ab -= (ps_aw * ps3_bw + ps_bw * ps3_aw + ps2_aw * ps2_bw) / ps_ww;
+            p3_ab += (ps_aw * ps2_bw * ps2_ww + ps_bw * ps2_aw * ps2_ww +
+                      ps_aw * ps_bw * ps3_ww) /
+              (ps_ww * ps_ww);
+          }
           gsl_matrix_set(PPPab, p, index_ab, p3_ab);
         }
       }
@@ -1187,21 +1192,23 @@ 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");
+
   gsl_vector *u_a = gsl_vector_safe_alloc(Uty->size);
 
-  for (size_t a = 1; a <= n_cvt + 2; ++a) {
+  for (size_t a = 1; a <= n_cvt + 2; ++a) { // walk columns of pheno+cvt
     if (a == n_cvt + 1) {
       continue;
     }
 
     if (a == n_cvt + 2) {
-      gsl_vector_safe_memcpy(u_a, Uty);
+      gsl_vector_safe_memcpy(u_a, Uty); // last column is phenotype
     } else {
       gsl_vector_const_view UtW_col = gsl_matrix_const_column(UtW, a - 1);
       gsl_vector_safe_memcpy(u_a, &UtW_col.vector);
     }
 
-    for (size_t b = a; b >= 1; --b) {
+    for (size_t b = a; b >= 1; --b) { // back fill other columns
       if (b == n_cvt + 1) {
         continue;
       }
@@ -1218,6 +1225,8 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
 
       gsl_vector_mul(&Uab_col.vector, u_a);
     }
+    cout << "a" << a << endl;
+    write(Uab,"Uab iteration");
   }
 
   gsl_vector_safe_free(u_a);
@@ -1963,11 +1972,16 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
   size_t n_cvt = UtW->size2, ni_test = UtW->size1;
   size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
 
+  cout << "n_cvt " << n_cvt << ", ni_test " << ni_test << ", n_index " << n_index << endl;
+
   gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index);
   gsl_vector *ab = gsl_vector_safe_alloc(n_index);
 
   gsl_matrix_set_zero(Uab);
+  write(UtW,"UtW");
+  write(UtW,"Uty");
   CalcUab(UtW, Uty, Uab);
+  write(Uab,"Uab");
   Calcab(UtW, Uty, ab);
 
   FUNC_PARAM param0 = {true, ni_test, n_cvt, eval, Uab, ab, 0};
diff --git a/src/param.cpp b/src/param.cpp
index 68e9d63..12f4299 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -1364,9 +1364,9 @@ void compAKtoS(const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt,
   * 1,1:0
   * 1,2:1
   * 1,3:2
-  * 2,2:3
-  * 2,3:4
-  * 3,3:5
+  * 2,2:5
+  * 2,3:6
+  * 3,3:9
 
 which is really the iteration moving forward along the diagonal and
 items to the right of it.
@@ -1385,8 +1385,8 @@ size_t GetabIndex(const size_t a, const size_t b, const size_t n_cvt) {
   size_t index = (2 * cols - a1 + 2) * (a1 - 1) / 2 + b1 - a1;
   cout << "* GetabIndx " << a1 << "," << b1 << "," << cols << ":" << index << endl;
   return index;
+  // return ( b < a ?  ((2 * cols - b + 2) * (b - 1) / 2 + a - b ): ((2 * cols - a + 2) * (a - 1) / 2 + b - a) );
 
-  // return ( b < a ?  ((2 * n - b + 2) * (b - 1) / 2 + a - b ): ((2 * n - a + 2) * (a - 1) / 2 + b - a) );
 }
 
 // From an existing n by nd (centered) G matrix, compute the d+1 by
generated by cgit-pink 1.4.1 (git 2.36.1) at 2025-09-16 12:27:19 +0000