Further debugging

6 files changed, 60 insertions, 53 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
diff --git a/test/src/unittests-math.cpp b/test/src/unittests-math.cpp
index 66a2f28..d44d8c4 100644
--- a/test/src/unittests-math.cpp
+++ b/test/src/unittests-math.cpp
@@ -12,34 +12,18 @@
 #include "mathfunc.h"
 #include "fastblas.h"
 #include "fastopenblas.h"
+#include "param.h"
 
 using namespace std;
 
 TEST_CASE("calcPab", "[math]") {
-// Hi_eval
-// vector size: 49
-  double Hi_eval[] = {1.0,0.00340397,0.00316875,0.00285475,0.000338947,0.000157865,0.000132286,9.72953e-05,8.30921e-05,7.85302e-05,7.50248e-05,6.77641e-05,6.74403e-05,5.97916e-05,3.66931e-05,3.56332e-05,3.26649e-05,2.94286e-05,2.78745e-05,2.57487e-05,2.34702e-05,2.23381e-05,2.17387e-05,2.11297e-05,1.99572e-05,1.92797e-05,1.79017e-05,1.77354e-05,1.76441e-05,1.69226e-05,1.6689e-05,1.62749e-05,1.56035e-05,1.53859e-05,1.50713e-05,1.46394e-05,1.43762e-05,1.35848e-05,1.29509e-05,9.6414e-06,9.07189e-06,8.36172e-06,7.36353e-06,7.28662e-06,6.93901e-06,4.85174e-06,3.93325e-06,3.44085e-06,1.68874e-06};
-// Uab
-// matrix size: 49 cols, 6 rows
-  double Uab[] = {49,8.88476e-27,1.60576e-26,2.79763e-26,1.04467e-27,8.13814e-28,3.43994e-29,1.83964e-31,1.00874e-29,1.51195e-29,5.16275e-29,2.30463e-30,3.49053e-30,6.02949e-30,1.02633e-30,4.2725e-29,1.29895e-30,2.51337e-31,2.0631e-30,1.35086e-32,5.46109e-31,6.47883e-31,1.80278e-30,2.23342e-30,1.89715e-30,4.93038e-30,8.54899e-31,1.67604e-31,2.68166e-30,1.11242e-30,4.93038e-32,2.6366e-31,3.15544e-30,2.41589e-30,7.64402e-31,4.77751e-32,2.49601e-31,1.92593e-32,2.94045e-30,1.70175e-30,1.17174e-30,5.117e-30,1.30193e-31,5.69767e-30,7.70372e-32,8.49335e-32,3.64434e-31,9.53456e-32,1.05464e-30,// row 0
-0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,// row 1
-0.0829848,-8.68273e-17,7.90197e-16,5.28651e-17,-3.03831e-16,3.1215e-17,-4.35439e-17,3.44317e-18,-2.91753e-17,-1.41231e-17,5.76483e-17,8.20483e-18,-1.90286e-18,1.4879e-17,-2.41261e-18,-1.20099e-16,-1.50572e-17,1.22551e-17,-1.98056e-17,1.29155e-19,-1.93584e-17,-1.16457e-17,-1.01071e-17,4.18941e-18,-2.53464e-17,2.53682e-17,3.17125e-18,-1.39346e-17,3.12915e-17,2.82287e-18,-4.43288e-18,-2.77284e-18,6.938e-18,2.71779e-17,2.64918e-17,-4.61261e-18,-5.52614e-18,1.17298e-18,-6.61135e-18,6.76754e-17,-1.20832e-17,-4.93118e-17,-4.31897e-18,-2.83118e-17,-3.32353e-18,2.9943e-19,-6.21337e-18,-1.92731e-17,4.96634e-17,// row 2
-0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,// row 3
-0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,// row 4
-0.00014054,8.48529e-07,3.88858e-05,9.98962e-08,8.83656e-05,1.19729e-06,5.51193e-05,6.44442e-05,8.43825e-05,1.31924e-05,6.43714e-05,2.92104e-05,1.03734e-06,3.67172e-05,5.67139e-06,0.000337593,0.00017454,0.000597554,0.000190132,1.23484e-06,0.000686213,0.000209331,5.66639e-05,7.85841e-06,0.000338635,0.000130526,1.17637e-05,0.00115853,0.000365131,7.1633e-06,0.000398558,2.91612e-05,1.52549e-05,0.000305741,0.000918125,0.00044534,0.000122349,7.14401e-05,1.4865e-05,0.00269132,0.000124604,0.000475211,0.000143276,0.000140681,0.000143383,1.05563e-06,0.000105934,0.00389584,0.00233867,// row 5
-  };
-// ab
-// vector size: 6
-  double ab[] = {6.91831e-310,6.91831e-310,0,0,0,0};
-// Pab
-// matrix size: 3 cols, 6 rows
-  double Pab[] = {49,4.08372e-33,1.2538e-76,// row 0
-0,1.72345e-47,4.65925e-33,// row 1
-0.0829848,3.06246e-57,3.2526e-86,// row 2
-0,0,1.36901e-71,// row 3
-0,0,4.17633e-62,// row 4
-0.000140915,-0.00014054,-0.00014054,// row 5
-  };
+  REQUIRE(GetabIndex(1,1,3) == 0);
+  REQUIRE(GetabIndex(1,2,3) == 1);
+  REQUIRE(GetabIndex(1,3,3) == 2);
+  REQUIRE(GetabIndex(2,2,3) == 5);
+  REQUIRE(GetabIndex(2,3,3) == 6);
+  REQUIRE(GetabIndex(3,3,3) == 9);
+  REQUIRE(GetabIndex(2,1,3) == 1); // reverse
 }
 
 TEST_CASE( "Math functions", "[math]" ) {