diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/debug.cpp | 19 | ||||
| -rw-r--r-- | src/debug.h | 5 | ||||
| -rw-r--r-- | src/lapack.cpp | 2 | ||||
| -rw-r--r-- | src/lmm.cpp | 19 | ||||
| -rw-r--r-- | src/param.cpp | 11 |
5 files changed, 44 insertions, 12 deletions
diff --git a/src/debug.cpp b/src/debug.cpp index 7728d92..37e4e8a 100644 --- a/src/debug.cpp +++ b/src/debug.cpp @@ -142,6 +142,7 @@ void disable_segfpe() { if (is_legacy_mode()) return; #ifdef __GNUC__ #if defined(__x86_64__) + debug_msg("disable segfpe"); fedisableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW); #endif #endif @@ -228,6 +229,24 @@ gsl_vector *gsl_vector_safe_alloc(size_t n) { return v; } +double do_gsl_matrix_safe_get(const gsl_matrix * m, const size_t i, const size_t j, + const char *__pretty_function, const char *__file, int __line) { + enforce(m); + if (!is_legacy_mode() && (is_debug_mode() || is_check_mode() || is_strict_mode())) { + auto rows = m->size1; + auto cols = m->size2; + if (i >= cols || j >= rows) { + std::string msg = "Matrix out of bounds (" + std::to_string(cols) + "," + std::to_string(rows) + ") "; + msg += std::to_string(i); + msg += ","; + msg += std::to_string(j); + fail_at_msg(__file,__line,msg.c_str()); + } + } + return gsl_matrix_get(m,i,j); +} + + char *do_strtok_safe(char *tokenize, const char *delimiters, const char *__pretty_function, const char *__file, int __line, const char *infile) { auto token = strtok(tokenize,delimiters); diff --git a/src/debug.h b/src/debug.h index b08d75d..64ab721 100644 --- a/src/debug.h +++ b/src/debug.h @@ -57,6 +57,9 @@ int gsl_matrix_safe_memcpy (gsl_matrix *dest, const gsl_matrix *src); void gsl_matrix_safe_free (gsl_matrix *v); void do_gsl_matrix_safe_free (gsl_matrix *m, const char *__pretty_function, const char *__file, int __line); +double do_gsl_matrix_safe_get (const gsl_matrix * m, const size_t i, const size_t j, const char *__pretty_function, const char *__file, int __line); +#define gsl_matrix_safe_get(m,i,j) do_gsl_matrix_safe_get(m, i, j,__SHOW_FUNC,__FILE__,__LINE__); + gsl_vector *gsl_vector_safe_alloc(size_t n); int gsl_vector_safe_memcpy (gsl_vector *dest, const gsl_vector *src); void gsl_vector_safe_free (gsl_vector *v); @@ -85,7 +88,7 @@ inline void warnfail_at_msg(bool strict, const char *__function, const char *__f } inline void fail_at_msg(const char *__file, int __line, std::string msg) { - std::cerr << msg << " in " << __file << " at line " << __line << std::endl; + std::cerr << "**** FAILED: " << msg << " in " << __file << " at line " << __line << std::endl; exit(1); } diff --git a/src/lapack.cpp b/src/lapack.cpp index 0fc2de4..6e43bd9 100644 --- a/src/lapack.cpp +++ b/src/lapack.cpp @@ -222,7 +222,7 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, if (INFO != 0) cerr << "ERROR: value of INFO is " << INFO; enforce_msg(INFO == 0, "lapack_eigen_symmv failed"); - if (is_check_mode()) disable_segfpe(); // reinstate fast NaN checking + if (is_check_mode()) enable_segfpe(); // reinstate fast NaN checking gsl_matrix_transpose(evec); diff --git a/src/lmm.cpp b/src/lmm.cpp index ed1366d..cdddf40 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -41,8 +41,6 @@ #include "gsl/gsl_roots.h" #include "gsl/gsl_vector.h" -// #include "eigenlib.h" - #include "gzstream.h" #include "gemma_io.h" #include "fastblas.h" @@ -226,9 +224,9 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval, double p_ab; double ps_ab, ps_aw, ps_bw, ps_ww; - for (size_t p = 0; p <= n_cvt + 1; ++p) { - for (size_t a = p + 1; a <= n_cvt + 2; ++a) { - for (size_t b = a; b <= n_cvt + 2; ++b) { + for (size_t p = 0; p <= n_cvt + 1; ++p) { // p walks phenotypes + covariates + for (size_t a = p + 1; a <= n_cvt + 2; ++a) { // a in p+1..rest + for (size_t b = a; b <= n_cvt + 2; ++b) { // b in a..rest index_ab = GetabIndex(a, b, n_cvt); if (p == 0) { gsl_vector_const_view Uab_col = @@ -244,11 +242,12 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval, index_bw = GetabIndex(b, p, n_cvt); index_ww = GetabIndex(p, p, n_cvt); - ps_ab = gsl_matrix_get(Pab, p - 1, index_ab); - ps_aw = gsl_matrix_get(Pab, p - 1, index_aw); - ps_bw = gsl_matrix_get(Pab, p - 1, index_bw); - ps_ww = gsl_matrix_get(Pab, p - 1, index_ww); + ps_ab = gsl_matrix_safe_get(Pab, p - 1, index_ab); + ps_aw = gsl_matrix_safe_get(Pab, p - 1, index_aw); + ps_bw = gsl_matrix_safe_get(Pab, p - 1, index_bw); + ps_ww = gsl_matrix_safe_get(Pab, p - 1, index_ww); + assert(ps_ww != 0.0); p_ab = ps_ab - ps_aw * ps_bw / ps_ww; gsl_matrix_set(Pab, p, index_ab, p_ab); } @@ -291,6 +290,7 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode, ps2_bw = gsl_matrix_get(PPab, p - 1, index_bw); ps2_ww = gsl_matrix_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; gsl_matrix_set(PPab, p, index_ab, p2_ab); @@ -340,6 +340,7 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode, ps3_bw = gsl_matrix_get(PPPab, p - 1, index_bw); ps3_ww = gsl_matrix_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; diff --git a/src/param.cpp b/src/param.cpp index edee79d..1a2b57c 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -1341,6 +1341,8 @@ void compAKtoS(const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt, if (tr_A == 0 || tr_K == 0) { d = 0; } else { + assert((tr_A * tr_K) - 1 != 0); + assert(ni_test - n_cvt != 0); d = d / (tr_A * tr_K) - 1 / (double)(ni_test - n_cvt); } @@ -1369,9 +1371,12 @@ size_t GetabIndex(const size_t a, const size_t b, const size_t n_cvt) { } size_t n = n_cvt + 2; - index = (2 * n - l + 2) * (l - 1) / 2 + h - l; + assert(2+h-1 != 0); + index = (2 * n - l + 2) * (l - 1) / 2 + h - l; return index; + + // 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 @@ -1386,6 +1391,7 @@ void compKtoV(const gsl_matrix *G, gsl_matrix *V) { gsl_vector *trKiKj = gsl_vector_alloc(n_vc * (n_vc + 1) / 2); gsl_vector *trKi = gsl_vector_alloc(n_vc); + assert(ni_test != 1); double d, tr, r = (double)ni_test / (double)(ni_test - 1); size_t t, t_il, t_jm, t_lm, t_im, t_jl, t_ij; @@ -1541,6 +1547,7 @@ void compKtoV(const gsl_matrix *G, gsl_matrix *V) { } } + assert(ni_test != 0); gsl_matrix_scale(V, 1.0 / pow((double)ni_test, 2)); gsl_matrix_free(KiKj); @@ -1603,6 +1610,7 @@ void JackknifeAKtoS(const gsl_matrix *W, const gsl_matrix *A, } for (size_t t = 0; t < ni_test; t++) { + assert(ni_test != 1); sumA[i][t] /= (double)(ni_test - 1); sumK[i][t] /= (double)(ni_test - 1); } @@ -1636,6 +1644,7 @@ void JackknifeAKtoS(const gsl_matrix *W, const gsl_matrix *A, } for (size_t t = 0; t < ni_test; t++) { + assert(ni_test != 1); sumAK[i][j][t] /= (double)(ni_test - 1); } |