Merge pull request #68 from genenetwork/issue26.

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

1 files changed, 34 insertions, 3 deletions

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++) {