Tests and enforces added related to https://github.com/genetics-statistics/GEMMA/issues/78

4 files changed, 47 insertions, 52 deletions

diff --git a/src/debug.h b/src/debug.h
index f47cb4c..e910a25 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -36,8 +36,8 @@ inline void fail_at_msg(bool strict, const char *__file, int __line, const char
 
 #else // DEBUG
 
-#define warning_msg(msg) cerr << "**** WARNING: " << msg << " in " << __FILE__ << " at line " << __LINE__ << " in " << __PRETTY_FUNCTION__ << endl;
-#define debug_msg(msg) cerr << "**** DEBUG: " << msg << " in " << __FILE__ << " at line " << __LINE__ << " in " << __PRETTY_FUNCTION__ << endl;
+#define warning_msg(msg) cerr << "**** WARNING: " << msg << " in " << __FILE__ << " at line " << __LINE__ << " in " << __FUNCTION__ << endl;
+#define debug_msg(msg) cerr << "**** DEBUG: " << msg << " in " << __FILE__ << " at line " << __LINE__ << " in " << __FUNCTION__ << endl;
 #define assert_issue(is_issue, expr) \
   ((is_issue) ? enforce_msg(expr,"FAIL: ISSUE assert") : __ASSERT_VOID_CAST(0))
 
diff --git a/src/io.cpp b/src/io.cpp
index fedf69a..f5783cd 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -133,7 +133,7 @@ std::istream &safeGetline(std::istream &is, std::string &t) {
   }
 }
 
-// Read SNP file.
+// Read SNP file. A single column of SNP names.
 bool ReadFile_snps(const string file_snps, set<string> &setSnps) {
   setSnps.clear();
 
@@ -148,6 +148,7 @@ bool ReadFile_snps(const string file_snps, set<string> &setSnps) {
 
   while (getline(infile, line)) {
     ch_ptr = strtok((char *)line.c_str(), " , \t");
+    enforce_msg(ch_ptr,"Problem reading SNP file");
     setSnps.insert(ch_ptr);
   }
 
@@ -157,6 +158,9 @@ bool ReadFile_snps(const string file_snps, set<string> &setSnps) {
   return true;
 }
 
+// Read SNP file using a header. The header determines how the
+// values for each row are parsed. A valid header can be, for example,
+// RS POS CHR
 bool ReadFile_snps_header(const string &file_snps, set<string> &setSnps) {
   setSnps.clear();
 
@@ -175,7 +179,7 @@ bool ReadFile_snps_header(const string &file_snps, set<string> &setSnps) {
   ReadHeader_io(line, header);
 
   if (header.rs_col == 0 && (header.chr_col == 0 || header.pos_col == 0)) {
-    cout << "missing rs id in the hearder" << endl;
+    cout << "missing rs id in the header" << endl;
   }
 
   while (!safeGetline(infile, line).eof()) {
@@ -185,6 +189,7 @@ bool ReadFile_snps_header(const string &file_snps, set<string> &setSnps) {
     ch_ptr = strtok((char *)line.c_str(), " , \t");
 
     for (size_t i = 0; i < header.coln; i++) {
+      enforce_msg(ch_ptr,"Problem reading SNP file");
       if (header.rs_col != 0 && header.rs_col == i + 1) {
         rs = ch_ptr;
       }
@@ -250,7 +255,7 @@ bool ReadFile_log(const string &file_log, double &pheno_mean) {
   return true;
 }
 
-// Read bimbam annotation file.
+// Read bimbam annotation file which consists of rows of SNP, POS and CHR
 bool ReadFile_anno(const string &file_anno, map<string, string> &mapRS2chr,
                    map<string, long int> &mapRS2bp,
                    map<string, double> &mapRS2cM) {
@@ -264,33 +269,39 @@ bool ReadFile_anno(const string &file_anno, map<string, string> &mapRS2chr,
   }
 
   string line;
-  char *ch_ptr;
-
-  string rs;
-  long int b_pos;
-  string chr;
-  double cM;
 
   while (!safeGetline(infile, line).eof()) {
-    ch_ptr = strtok((char *)line.c_str(), " , \t");
-    rs = ch_ptr;
+    const char *ch_ptr = strtok((char *)line.c_str(), " , \t");
+    enforce_str(ch_ptr, line + " Bad RS format");
+    const string rs = ch_ptr;
+    enforce_str(rs != "", line + " Bad RS format");
+
     ch_ptr = strtok(NULL, " , \t");
+    enforce_str(ch_ptr, line + " Bad format");
+    ulong b_pos;
     if (strcmp(ch_ptr, "NA") == 0) {
       b_pos = -9;
     } else {
       b_pos = atol(ch_ptr);
     }
+    enforce_str(b_pos,line + " Bad pos format (is zero)");
+
+    string chr;
     ch_ptr = strtok(NULL, " , \t");
     if (ch_ptr == NULL || strcmp(ch_ptr, "NA") == 0) {
       chr = "-9";
     } else {
       chr = ch_ptr;
+      enforce_str(chr != "", line + " Bad chr format");
     }
+
+    double cM;
     ch_ptr = strtok(NULL, " , \t");
     if (ch_ptr == NULL || strcmp(ch_ptr, "NA") == 0) {
       cM = -9;
     } else {
       cM = atof(ch_ptr);
+      enforce_str(b_pos, line + "Bad cM format (is zero)");
     }
 
     mapRS2chr[rs] = chr;
@@ -377,9 +388,9 @@ bool ReadFile_pheno(const string &file_pheno,
 
   while (!safeGetline(infile, line).eof()) {
     ch_ptr = strtok((char *)line.c_str(), " , \t");
-
     size_t i = 0;
     while (i < p_max) {
+      enforce_msg(ch_ptr,"Wrong number of phenotypes");
       if (mapP2c.count(i + 1) != 0) {
         if (strcmp(ch_ptr, "NA") == 0) {
           ind_pheno_row[mapP2c[i + 1]] = 0;
@@ -618,7 +629,8 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
   gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
   int sig;
   LUDecomp(WtW, pmt, &sig);
-  LUInvert(WtW, pmt, WtWi);
+
+  LUInvert(WtW, pmt, WtWi); // @@
 
   double v_x, v_w;
   int c_idv = 0;
@@ -667,6 +679,10 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
     }
 
     if (mapRS2bp.count(rs) == 0) {
+      std::string msg = "Can't figure out position for ";
+      msg += rs;
+      debug_msg(msg);
+
       chr = "-9";
       b_pos = -9;
       cM = -9;
@@ -2979,7 +2995,7 @@ bool bgenKin(const string &file_oxford, vector<int> &indicator_snp,
 bool ReadHeader_io(const string &line, HEADER &header) {
   string rs_ptr[] = {"rs",    "RS",    "snp",  "SNP",  "snps",      "SNPS",
                      "snpid", "SNPID", "rsid", "RSID", "MarkerName"};
-  set<string> rs_set(rs_ptr, rs_ptr + 11);
+  set<string> rs_set(rs_ptr, rs_ptr + 11); // create a set of 11 items
   string chr_ptr[] = {"chr", "CHR"};
   set<string> chr_set(chr_ptr, chr_ptr + 2);
   string pos_ptr[] = {
diff --git a/src/lapack.cpp b/src/lapack.cpp
index b7c6f74..9e25791 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -24,6 +24,7 @@
 #include <vector>
 
 #include "debug.h"
+#include "lapack.h"
 #include "mathfunc.h"
 
 using namespace std;
@@ -520,7 +521,7 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty,
 
 // LU decomposition.
 void LUDecomp(gsl_matrix *LU, gsl_permutation *p, int *signum) {
-  gsl_linalg_LU_decomp(LU, p, signum);
+  enforce_gsl(gsl_linalg_LU_decomp(LU, p, signum));
   return;
 }
 
@@ -551,11 +552,13 @@ void LUDecomp(gsl_matrix_float *LU, gsl_permutation *p, int *signum) {
 }
 */
 
-// LU invert.
-void LUInvert(const gsl_matrix *LU, const gsl_permutation *p,
-              gsl_matrix *inverse) {
-  gsl_linalg_LU_invert(LU, p, inverse);
-  return;
+// LU invert. Returns inverse. Note that GSL does not recommend using
+// this function
+void LUInvert(const gsl_matrix *LU, const gsl_permutation *p, gsl_matrix *ret_inverse) {
+  auto det = LULndet(LU);
+  assert(det != 1.0);
+
+  enforce_gsl(gsl_linalg_LU_invert(LU, p, ret_inverse));
 }
 
 /*
@@ -590,9 +593,10 @@ void LUInvert(const gsl_matrix_float *LU, const gsl_permutation *p,
 */
 
 // LU lndet.
-double LULndet(gsl_matrix *LU) {
+double LULndet(const gsl_matrix *LU) {
   double d;
-  d = gsl_linalg_LU_lndet(LU);
+  d = gsl_linalg_LU_lndet((gsl_matrix *)LU);
+  enforce(d != 0.0);
   return d;
 }
 
@@ -620,7 +624,7 @@ double LULndet(gsl_matrix_float *LU) {
 // LU solve.
 void LUSolve(const gsl_matrix *LU, const gsl_permutation *p,
              const gsl_vector *b, gsl_vector *x) {
-  gsl_linalg_LU_solve(LU, p, b, x);
+  enforce_gsl(gsl_linalg_LU_solve(LU, p, b, x));
   return;
 }
 
diff --git a/src/lapack.h b/src/lapack.h
index 9596667..9115201 100644
--- a/src/lapack.h
+++ b/src/lapack.h
@@ -23,47 +23,22 @@
 
 using namespace std;
 
-void lapack_float_cholesky_decomp(gsl_matrix_float *A);
 void lapack_cholesky_decomp(gsl_matrix *A);
-void lapack_float_cholesky_solve(gsl_matrix_float *A, const gsl_vector_float *b,
-                                 gsl_vector_float *x);
 void lapack_cholesky_solve(gsl_matrix *A, const gsl_vector *b, gsl_vector *x);
-void lapack_sgemm(char *TransA, char *TransB, float alpha,
-                  const gsl_matrix_float *A, const gsl_matrix_float *B,
-                  float beta, gsl_matrix_float *C);
 void lapack_dgemm(char *TransA, char *TransB, double alpha, const gsl_matrix *A,
                   const gsl_matrix *B, double beta, gsl_matrix *C);
-void lapack_float_eigen_symmv(gsl_matrix_float *A, gsl_vector_float *eval,
-                              gsl_matrix_float *evec,
-                              const size_t flag_largematrix);
 void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
                         const size_t flag_largematrix);
-
 double EigenDecomp(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
                    const size_t flag_largematrix);
 double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
                    const size_t flag_largematrix);
-double EigenDecomp(gsl_matrix_float *G, gsl_matrix_float *U,
-                   gsl_vector_float *eval, const size_t flag_largematrix);
-
 double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty);
-double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty,
-                     gsl_vector_float *OiXty);
-
 void LUDecomp(gsl_matrix *LU, gsl_permutation *p, int *signum);
-void LUDecomp(gsl_matrix_float *LU, gsl_permutation *p, int *signum);
-void LUInvert(const gsl_matrix *LU, const gsl_permutation *p,
-              gsl_matrix *inverse);
-void LUInvert(const gsl_matrix_float *LU, const gsl_permutation *p,
-              gsl_matrix_float *inverse);
-double LULndet(gsl_matrix *LU);
-double LULndet(gsl_matrix_float *LU);
+void LUInvert(const gsl_matrix *LU, const gsl_permutation *p, gsl_matrix *ret_inverse);
+double LULndet(const gsl_matrix *LU);
 void LUSolve(const gsl_matrix *LU, const gsl_permutation *p,
              const gsl_vector *b, gsl_vector *x);
-void LUSolve(const gsl_matrix_float *LU, const gsl_permutation *p,
-             const gsl_vector_float *b, gsl_vector_float *x);
-
 bool lapack_ddot(vector<double> &x, vector<double> &y, double &v);
-bool lapack_sdot(vector<float> &x, vector<float> &y, double &v);
 
 #endif