diff options
| -rw-r--r-- | src/debug.cpp | 35 | ||||
| -rw-r--r-- | src/debug.h | 3 | ||||
| -rw-r--r-- | src/gemma.cpp | 17 | ||||
| -rwxr-xr-x | test/dev_test_suite.sh | 5 | ||||
| -rwxr-xr-x | test/test_suite.sh | 12 |
5 files changed, 65 insertions, 7 deletions
diff --git a/src/debug.cpp b/src/debug.cpp index 4e58d5d..0b9c19d 100644 --- a/src/debug.cpp +++ b/src/debug.cpp @@ -26,6 +26,7 @@ #include <string> #include <sys/stat.h> #include <vector> +#include <fenv.h> #include "gsl/gsl_blas.h" #include "gsl/gsl_cdf.h" @@ -59,6 +60,40 @@ bool is_quiet_mode() { return debug_quiet; }; bool is_issue(uint issue) { return issue == debug_issue; }; bool is_legacy_mode() { return debug_legacy; }; +#include <stdio.h> +#include <sys/types.h> +#include <unistd.h> +#include <signal.h> + +void sighandler(int signum) +{ + cout << R"( +FATAL ERROR: GEMMA caused a floating point error which suggests machine boundaries were reached. + +You can disable floating point tests with the -no-check switch (use at your own risk!) +)" << endl; + signal(signum, SIG_DFL); + kill(getpid(), signum); // should force a core dump +} + +/* + Force the floating point processor to throw an exception when the result of + a double/float computation is overflow, underflow, NaN or inf. In principle + this is an Intel hardware feature that does not slow down computations. +*/ + +void enable_segfpe() { + #ifdef __GNUC__ + signal(SIGFPE, sighandler); + feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW); + #endif +} + +void disable_segfpe() { + #ifdef __GNUC__ + fedisableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW); + #endif +} /* Helper function to make sure gsl allocations do their job because diff --git a/src/debug.h b/src/debug.h index 67764df..ced8e72 100644 --- a/src/debug.h +++ b/src/debug.h @@ -45,6 +45,9 @@ bool is_quiet_mode(); bool is_issue(uint issue); bool is_legacy_mode(); +void enable_segfpe(); +void disable_segfpe(); + #define check_int_mult_overflow(m,n) \ { auto x = m * n; \ enforce_msg(x / m == n, "multiply integer overflow"); } diff --git a/src/gemma.cpp b/src/gemma.cpp index 809cd8e..a183edc 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -1624,6 +1624,8 @@ void GEMMA::BatchRun(PARAM &cPar) { clock_t time_begin, time_start; time_begin = clock(); + if (is_check_mode()) enable_segfpe(); // fast NaN checking + // Read Files. cout << "Reading Files ... " << endl; cPar.ReadFiles(); @@ -1882,6 +1884,7 @@ void GEMMA::BatchRun(PARAM &cPar) { enforce_msg(G, "allocate G"); // just to be sure time_start = clock(); + cPar.CalcKin(G); cPar.time_G = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0); @@ -1909,6 +1912,8 @@ void GEMMA::BatchRun(PARAM &cPar) { VARCOV cVarcov; cVarcov.CopyFromParam(cPar); + if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (!cPar.file_bfile.empty()) { cVarcov.AnalyzePlink(); } else { @@ -2022,6 +2027,8 @@ void GEMMA::BatchRun(PARAM &cPar) { VARCOV cVarcov; cVarcov.CopyFromParam(cPar); + if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (!cPar.file_bfile.empty()) { cVarcov.AnalyzePlink(); } else { @@ -2048,6 +2055,8 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_vector_view Y_col = gsl_matrix_column(Y, 0); + if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (!cPar.file_gene.empty()) { cLm.AnalyzeGene(W, &Y_col.vector); // y is the predictor, not the phenotype @@ -2754,6 +2763,8 @@ void GEMMA::BatchRun(PARAM &cPar) { LMM cLmm; cLmm.CopyFromParam(cPar); + if (is_check_mode()) disable_segfpe(); // fast NaN checking + gsl_vector_view Y_col = gsl_matrix_column(Y, 0); gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); @@ -2770,6 +2781,7 @@ void GEMMA::BatchRun(PARAM &cPar) { } else { // BIMBAM analysis + if (cPar.file_gxe.empty()) { cLmm.AnalyzeBimbam(U, eval, UtW, &UtY_col.vector, W, &Y_col.vector, cPar.setGWASnps); @@ -2785,6 +2797,8 @@ void GEMMA::BatchRun(PARAM &cPar) { MVLMM cMvlmm; cMvlmm.CopyFromParam(cPar); + if (is_check_mode()) disable_segfpe(); // fast NaN checking + if (!cPar.file_bfile.empty()) { if (cPar.file_gxe.empty()) { cMvlmm.AnalyzePlink(U, eval, UtW, UtY); @@ -3107,6 +3121,9 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) { outfile << "##" << endl; outfile << "## GEMMA Version = " << version << " (" << date << ")" << endl; outfile << "## Build profile = " << GEMMA_PROFILE << endl ; + #ifdef __GNUC__ + outfile << "## GCC version = " << __GNUC__ << "." << __GNUC_MINOR__ << "." << __GNUC_PATCHLEVEL__ << endl; + #endif outfile << "## GSL Version = " << GSL_VERSION << endl; outfile << "## Eigen Version = " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << endl; #ifdef OPENBLAS diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh index 0d3d8a0..ad743ff 100755 --- a/test/dev_test_suite.sh +++ b/test/dev_test_suite.sh @@ -14,6 +14,7 @@ testBXDStandardRelatednessMatrixKSingularError() { -c ../example/BXD_covariates.txt \ -a ../example/BXD_snps.txt \ -gk \ + -no-check \ -o $outn assertEquals 22 $? # should show singular error } @@ -56,7 +57,7 @@ testBXDLMMLikelihoodRatio() { -c ../example/BXD_covariates2.txt \ -a ../example/BXD_snps.txt \ -k ./output/BXD.cXX.txt \ - -lmm 2 -maf 0.1 \ + -lmm 2 -no-check -maf 0.1 \ -o $outn assertEquals 0 $? @@ -89,7 +90,7 @@ testUnivariateLinearMixedModelLOCO1() { -a ../example/mouse_hs1940.anno.txt \ -k ./output/mouse_hs1940_LOCO1.cXX.txt \ -snps ../example/mouse_hs1940_snps.txt -lmm \ - -nind 400 \ + -nind 400 -no-check \ -o $outn assertEquals 0 $? grep "total computation time" < output/$outn.log.txt diff --git a/test/test_suite.sh b/test/test_suite.sh index 7af33aa..faddd43 100755 --- a/test/test_suite.sh +++ b/test/test_suite.sh @@ -8,7 +8,7 @@ testBslmm1() { $gemma $gemmaopts -g ../example/mouse_hs1940.geno.txt.gz \ -p ../example/mouse_hs1940.pheno.txt \ -n 2 -a ../example/mouse_hs1940.anno.txt \ - -bslmm \ + -bslmm -no-check \ -o $outn -w 1000 -s 10000 -seed 1 assertEquals 0 $? outfn1=output/$outn.hyp.txt @@ -39,7 +39,7 @@ testBslmm3() { -a ../example/mouse_hs1940.anno.txt \ -bslmm \ -o $outn \ - -w 1000 -s 10000 -seed 1 + -w 1000 -s 10000 -seed 1 -no-check assertEquals 0 $? outfn=output/$outn.hyp.txt # assertEquals "291" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.0f",(substr($x,,0,6))) } END { printf "%.0f",100*$sum }' $outfn` @@ -54,7 +54,7 @@ testBslmm4() { -emu ./output/mouse_hs1940_CD8_bslmm.log.txt \ -ebv ./output/mouse_hs1940_CD8_bslmm.bv.txt \ -k ./output/mouse_hs1940_CD8_train.cXX.txt \ - -predict \ + -predict -no-check \ -o $outn assertEquals 0 $? outfn=output/$outn.prdt.txt @@ -98,7 +98,7 @@ testUnivariateLinearMixedModelFullLOCO1() { -loco 1 \ -a ../example/mouse_hs1940.anno.txt \ -k ./output/mouse_hs1940_full_LOCO1.cXX.txt \ - -lmm \ + -lmm -no-check \ -o $outn assertEquals 0 $? grep "total computation time" < output/$outn.log.txt @@ -142,7 +142,8 @@ testLinearMixedModelPhenotypes() { -n 1 6 \ -a ../example/mouse_hs1940.anno.txt \ -k ./output/mouse_hs1940.cXX.txt \ - -lmm -o mouse_hs1940_CD8MCH_lmm + -lmm -no-check \ + -o mouse_hs1940_CD8MCH_lmm assertEquals 0 $? outfn=output/mouse_hs1940_CD8MCH_lmm.assoc.txt @@ -172,6 +173,7 @@ testPlinkLinearMixedModelCovariates() { -lmm 1 \ -maf 0.1 \ -c $datadir/HLC_covariates.txt \ + -no-check \ -o $testname assertEquals 0 $? outfn=output/$testname.assoc.txt |