aboutsummaryrefslogtreecommitdiff
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 v1.2.3 (git 2.44.0) at 2025-01-18 19:27:46 +0000