Tests and fixes https://github.com/genetics-statistics/GEMMA/issues/26

6 files changed, 60 insertions, 3 deletions

diff --git a/src/debug.h b/src/debug.h
index fafebcf..3fbe9e0 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -7,10 +7,15 @@
 // enforce works like assert but also when NDEBUG is set (i.e., it
 // always works). enforce_msg prints message instead of expr
 
+#define ROUND(f) round(f * 10000.)/10000
 #if defined NDEBUG
 #define debug_msg(msg)
+#define assert_issue(is_issue, expr)
 #else
 #define debug_msg(msg) cout << "**** DEBUG: " << msg << endl;
+#define assert_issue(is_issue, expr) \
+  ((is_issue) ? enforce_msg(expr,"FAIL: ISSUE assert") : __ASSERT_VOID_CAST(0))
+
 #endif
 
 /* This prints an "Assertion failed" message and aborts.  */
@@ -49,4 +54,5 @@ inline void __enforce_fail(const char *__assertion, const char *__file,
        : __enforce_fail(gsl_strerror(COMBINE(res, __LINE__)), __FILE__,         \
                         __LINE__, __ASSERT_FUNCTION))
 
+
 #endif
diff --git a/src/gemma.cpp b/src/gemma.cpp
index a678b8c..3197043 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -697,6 +697,7 @@ void GEMMA::PrintHelp(size_t option) {
     cout << " -silence                 silent terminal display" << endl;
     cout << " -debug                   debug output" << endl;
     cout << " -nind       [num]        read up to num individuals" << endl;
+    cout << " -issue      [num]        enable tests relevant to issue tracker" << endl;
     cout << endl;
   }
 
@@ -1347,6 +1348,16 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) {
       str.clear();
       str.assign(argv[i]);
       cPar.ni_max = atoi(str.c_str()); // for testing purposes
+      enforce(cPar.ni_max > 0);
+    } else if (strcmp(argv[i], "-issue") == 0) {
+      if (argv[i + 1] == NULL || argv[i + 1][0] == '-') {
+        continue;
+      }
+      ++i;
+      str.clear();
+      str.assign(argv[i]);
+      cPar.issue = atoi(str.c_str()); // for testing purposes
+      enforce(cPar.issue > 0);
     } else if (strcmp(argv[i], "-emp") == 0) {
       if (argv[i + 1] == NULL || argv[i + 1][0] == '-') {
         continue;
@@ -2656,11 +2667,12 @@ return;}
     gsl_matrix *se_B = gsl_matrix_alloc(Y->size2, W->size2);
     gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1);
     gsl_matrix *U = gsl_matrix_alloc(Y->size1, Y->size1);
-    gsl_matrix *UtW = gsl_matrix_alloc(Y->size1, W->size2);
-    gsl_matrix *UtY = gsl_matrix_alloc(Y->size1, Y->size2);
-    gsl_vector *eval = gsl_vector_alloc(Y->size1);
+    gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2);
+    gsl_matrix *UtY = gsl_matrix_calloc(Y->size1, Y->size2);
+    gsl_vector *eval = gsl_vector_calloc(Y->size1);
     gsl_vector *env = gsl_vector_alloc(Y->size1);
     gsl_vector *weight = gsl_vector_alloc(Y->size1);
+    assert_issue(cPar.issue == 26, UtY->data[0] == 0.0);
 
     // set covariates matrix W and phenotype matrix Y
     // an intercept should be included in W,
@@ -2769,6 +2781,8 @@ return;}
       CalcUtX(U, W, UtW);
       CalcUtX(U, Y, UtY);
 
+      assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
+
       LMM cLmm;
       cLmm.CopyFromParam(cPar);
 
@@ -2784,6 +2798,7 @@ return;}
       // calculate UtW and Uty
       CalcUtX(U, W, UtW);
       CalcUtX(U, Y, UtY);
+      assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
 
       // calculate REMLE/MLE estimate and pve for univariate model
       if (cPar.n_ph == 1) { // one phenotype
@@ -2791,24 +2806,40 @@ return;}
         gsl_vector_view se_beta = gsl_matrix_row(se_B, 0);
         gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
 
+        assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143);
+
         CalcLambda('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max,
                    cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
+        assert(!std::isnan(UtY->data[0]));
+        assert(!std::isnan(B->data[0]));
+        assert(!std::isnan(se_B->data[0]));
+
         CalcLmmVgVeBeta(eval, UtW, &UtY_col.vector, cPar.l_mle_null,
                         cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector,
                         &se_beta.vector);
 
+        assert(!std::isnan(UtY->data[0]));
+        assert(!std::isnan(B->data[0]));
+        assert(!std::isnan(se_B->data[0]));
+
         cPar.beta_mle_null.clear();
         cPar.se_beta_mle_null.clear();
         for (size_t i = 0; i < B->size2; i++) {
           cPar.beta_mle_null.push_back(gsl_matrix_get(B, 0, i));
           cPar.se_beta_mle_null.push_back(gsl_matrix_get(se_B, 0, i));
         }
+        assert(!std::isnan(UtY->data[0]));
+        assert(!std::isnan(B->data[0]));
+        assert(!std::isnan(se_B->data[0]));
+        assert(!std::isnan(cPar.beta_mle_null.front()));
+        assert(!std::isnan(cPar.se_beta_mle_null.front()));
 
         CalcLambda('R', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max,
                    cPar.n_region, cPar.l_remle_null, cPar.logl_remle_H0);
         CalcLmmVgVeBeta(eval, UtW, &UtY_col.vector, cPar.l_remle_null,
                         cPar.vg_remle_null, cPar.ve_remle_null, &beta.vector,
                         &se_beta.vector);
+
         cPar.beta_remle_null.clear();
         cPar.se_beta_remle_null.clear();
         for (size_t i = 0; i < B->size2; i++) {
diff --git a/src/main.cpp b/src/main.cpp
index 833136c..e37b154 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -25,10 +25,20 @@
 
 using namespace std;
 
+void gemma_gsl_error_handler (const char * reason,
+                              const char * file,
+                              int line, int gsl_errno) {
+  cerr << "GSL ERROR: " << reason << " in " << file
+       << " at line " << line << " errno " << gsl_errno <<endl;
+  exit(22);
+}
+
 int main(int argc, char *argv[]) {
   GEMMA cGemma;
   PARAM cPar;
 
+  gsl_set_error_handler (&gemma_gsl_error_handler);
+
   if (argc <= 1) {
     cGemma.PrintHeader();
     return EXIT_SUCCESS;
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index 9e19bf1..9a4bb8b 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -187,6 +187,14 @@ double ScaleMatrix(gsl_matrix *G) {
   return d;
 }
 
+double SumVector(const gsl_vector *v) {
+  double sum = 0;
+  for (int i = 0; i < v->size; i++ ) {
+    sum += gsl_vector_get(v, i);
+  }
+  return( sum );
+}
+
 // Center the vector y.
 double CenterVector(gsl_vector *y) {
   double d = 0.0;
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 29eafe4..b9f3c73 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -32,6 +32,7 @@ void CenterMatrix(gsl_matrix *G, const gsl_vector *w);
 void CenterMatrix(gsl_matrix *G, const gsl_matrix *W);
 void StandardizeMatrix(gsl_matrix *G);
 double ScaleMatrix(gsl_matrix *G);
+double SumVector(const gsl_vector *v);
 double CenterVector(gsl_vector *y);
 void CenterVector(gsl_vector *y, const gsl_matrix *W);
 void StandardizeVector(gsl_vector *y);
diff --git a/src/param.h b/src/param.h
index 45d8c0f..f4d649f 100644
--- a/src/param.h
+++ b/src/param.h
@@ -114,6 +114,7 @@ public:
   // IO-related parameters.
   bool mode_silence;
   bool mode_debug = false;
+  uint issue; // enable tests for issue on github tracker
   int a_mode; // Analysis mode, 1/2/3/4 for Frequentist tests
   int k_mode; // Kinship read mode: 1: n by n matrix, 2: id/id/k_value;
   vector<size_t> p_column; // Which phenotype column needs analysis.