Add floating point hardware checking for Intel on GNU compilers

When using the -check function (the default) it is enabled for Kinship
computation and LM/LMM up to individual SNP computation. This means
there can no longer be NaN values for matrices that are reused for every
SNP, but it is possible to have NaN for individual SNPs.

Fixes #161

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

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"); }

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

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

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