segfpe handling

5 files changed, 30 insertions, 13 deletions

diff --git a/INSTALL.md b/INSTALL.md
index f4cc046..be4469f 100644
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -63,9 +63,10 @@ issues (i.e., for the purpose of debugging).
 
 Note that this is an advanced configuration option at this stage. GNU
 Guix will make it easier in the future to deal with shared
-graphs. Contact Pjotr Prins if you are really interested.
+graphs. Contact Pjotr Prins if you are really interested. For development
+we have the [gemma-dev-env package](https://gitlab.com/genenetwork/guix-bioinformatics/blob/master/gn/packages/gemma.scm).
 
-The following two links capture and provide the reproducible build
+The following two links capture and provide the current reproducible build
 system that we use for development of GEMMA to generate /gnu/store/9ahrb1swr06kjm2gr2zg0fsyvps3xqgz-profile
 
 1. https://gitlab.com/genenetwork/guix-bioinformatics/tree/99718d253ec9ed8ed836f0a348381a7cd83d4b9f
diff --git a/src/gemma.cpp b/src/gemma.cpp
index a183edc..78633e8 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1624,7 +1624,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
   clock_t time_begin, time_start;
   time_begin = clock();
 
-  if (is_check_mode()) enable_segfpe(); // fast NaN checking
+  if (is_check_mode()) enable_segfpe(); // fast NaN checking by default
 
   // Read Files.
   cout << "Reading Files ... " << endl;
@@ -1912,7 +1912,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     VARCOV cVarcov;
     cVarcov.CopyFromParam(cPar);
 
-    if (is_check_mode()) disable_segfpe(); // fast NaN checking
+    if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
     if (!cPar.file_bfile.empty()) {
       cVarcov.AnalyzePlink();
@@ -2027,7 +2027,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     VARCOV cVarcov;
     cVarcov.CopyFromParam(cPar);
 
-    if (is_check_mode()) disable_segfpe(); // fast NaN checking
+    if (is_check_mode()) disable_segfpe(); // fast NaN checking for now
 
     if (!cPar.file_bfile.empty()) {
       cVarcov.AnalyzePlink();
@@ -2055,7 +2055,7 @@ 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 (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
       if (!cPar.file_gene.empty()) {
         cLm.AnalyzeGene(W,
@@ -2763,7 +2763,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
           LMM cLmm;
           cLmm.CopyFromParam(cPar);
 
-          if (is_check_mode()) disable_segfpe(); // fast NaN checking
+          if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
           gsl_vector_view Y_col = gsl_matrix_column(Y, 0);
           gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
@@ -2790,14 +2790,13 @@ void GEMMA::BatchRun(PARAM &cPar) {
                                     &Y_col.vector, env);
             }
           }
-
           cLmm.WriteFiles();
           cLmm.CopyToParam(cPar);
         } else {
           MVLMM cMvlmm;
           cMvlmm.CopyFromParam(cPar);
 
-          if (is_check_mode()) disable_segfpe(); // fast NaN checking
+          if (is_check_mode()) disable_segfpe(); // disable fast NaN checking
 
           if (!cPar.file_bfile.empty()) {
             if (cPar.file_gxe.empty()) {
diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp
index 818c5e8..fad3e47 100644
--- a/src/gemma_io.cpp
+++ b/src/gemma_io.cpp
@@ -1483,12 +1483,14 @@ bool BimbamKin(const string file_geno, const set<string> ksnps,
   // and transpose
   // FIXME: the following is very slow
 
+  debug_msg("begin transpose");
   for (size_t i = 0; i < ni_total; ++i) {
     for (size_t j = 0; j < i; ++j) {
       d = gsl_matrix_get(matrix_kin, j, i);
       gsl_matrix_set(matrix_kin, i, j, d);
     }
   }
+  debug_msg("end transpose");
   // GSL is faster - and there are even faster methods
   // enforce_gsl(gsl_matrix_transpose(matrix_kin));
 
diff --git a/src/lapack.cpp b/src/lapack.cpp
index 187daa8..b309a1b 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -193,6 +193,8 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
     double WORK_temp[1];
     int IWORK_temp[1];
 
+    if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
+
     // DSYEVR - computes selected eigenvalues and, optionally,
     // eigenvectors of a real symmetric matrix
     // Here compute both (JOBZ is V), all eigenvalues (RANGE is A)
@@ -200,7 +202,12 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
     dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU,
             &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp,
             &LWORK, IWORK_temp, &LIWORK, &INFO);
-    enforce_msg(INFO != 0, "lapack_eigen_symmv failed");
+    // If info = 0, the execution is successful.
+    // If info = -i, the i-th parameter had an illegal value.
+    // If info = i, an internal error has occurred.
+
+    if (INFO != 0) cerr << "ERROR: value of INFO is " << INFO;
+    enforce_msg(INFO == 0, "lapack_eigen_symmv failed");
     LWORK = (int)WORK_temp[0];
     LIWORK = (int)IWORK_temp[0];
 
@@ -210,7 +217,10 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
     dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU,
             &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK,
             IWORK, &LIWORK, &INFO);
-    enforce_msg(INFO != 0, "lapack_eigen_symmv failed");
+    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
 
     gsl_matrix_transpose(evec);
 
diff --git a/src/lmm.cpp b/src/lmm.cpp
index acd9667..7e4c3e1 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -97,8 +97,9 @@ void LMM::CopyToParam(PARAM &cPar) {
 }
 
 void LMM::WriteFiles() {
-
   string file_str;
+  debug_msg("LMM::WriteFiles");
+
   file_str = path_out + "/" + file_out;
   file_str += ".assoc.txt";
 
@@ -1319,6 +1320,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
 
   auto batch_compute = [&](size_t l) { // using a C++ closure
     // Compute SNPs in batch, note the computations are independent per SNP
+    debug_msg("enter batch_compute");
     gsl_matrix_view Xlarge_sub = gsl_matrix_submatrix(Xlarge, 0, 0, inds, l);
     gsl_matrix_view UtXlarge_sub =
         gsl_matrix_submatrix(UtXlarge, 0, 0, inds, l);
@@ -1366,6 +1368,7 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
                       p_wald, p_lrt, p_score, logl_H1};
       sumStat.push_back(SNPs);
     }
+    debug_msg("exit batch_compute");
   };
 
   const auto num_snps = indicator_snp.size();
@@ -1441,9 +1444,11 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
     gsl_vector_safe_memcpy(&Xlarge_col.vector, x);
     c++; // count SNPs going in
 
-    if (c % msize == 0)
+    if (c % msize == 0) {
       batch_compute(msize);
+    }
   }
+
   batch_compute(c % msize);
   ProgressBar("Reading SNPs", num_snps - 1, num_snps - 1);
   // cout << "Counted SNPs " << c << " sumStat " << sumStat.size() << endl;