Using M_LLM? modes

3 files changed, 30 insertions, 21 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index 1938e56..5a6d4d8 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2550,8 +2550,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
   // LMM or mvLMM or Eigen-Decomposition
   if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 ||
-      cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 ||
-      cPar.a_mode == M_EIGEN) { // Fit LMM or mvLMM or eigen
+      cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 || cPar.a_mode == M_LMM9 ||
+      cPar.a_mode == M_EIGEN  ) { // Fit LMM or mvLMM or eigen
     gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph);
     enforce_msg(Y, "allocate Y"); // just to be sure
     gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);
@@ -2565,6 +2565,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     gsl_vector *eval = gsl_vector_calloc(Y->size1);
     gsl_vector *env = gsl_vector_safe_alloc(Y->size1);
     gsl_vector *weight = gsl_vector_safe_alloc(Y->size1);
+    debug_msg("Started on LMM");
     assert_issue(is_issue(26), UtY->data[0] == 0.0);
 
     // set covariates matrix W and phenotype matrix Y
@@ -2578,6 +2579,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     if (!(cPar.file_kin).empty()) {
       ReadFile_kin(cPar.file_kin, cPar.indicator_idv, cPar.mapID2num,
                    cPar.k_mode, cPar.error, G);
+      debug_msg("Read K/GRM file");
       if (cPar.error == true) {
         cout << "error! fail to read kinship/relatedness file. " << endl;
         return;
@@ -2609,7 +2611,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
         }
       }
 
-      // eigen-decomposition and calculate trace_G
+      // eigen-decomposition and calculate trace_G - main track
       cout << "Start Eigen-Decomposition..." << endl;
       time_start = clock();
 
@@ -2663,7 +2665,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     if (cPar.a_mode == M_EIGEN) {
       cPar.WriteMatrix(U, "eigenU");
       cPar.WriteVector(eval, "eigenD");
-    } else if (!cPar.file_gene.empty()) {
+    } else if (!cPar.file_gene.empty()) { // Run with gene file
       // calculate UtW and Uty
       CalcUtX(U, W, UtW);
       CalcUtX(U, Y, UtY);
@@ -2682,6 +2684,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cLmm.WriteFiles();
       cLmm.CopyToParam(cPar);
     } else {
+      debug_msg("Main LMM track");
       // calculate UtW and Uty
       CalcUtX(U, W, UtW);
       CalcUtX(U, Y, UtY);
@@ -2736,6 +2739,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
         CalcPve(eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.trace_G,
                 cPar.pve_null, cPar.pve_se_null);
+        debug_msg("main print summary");
         cPar.PrintSummary();
 
         // calculate and output residuals
@@ -2771,15 +2775,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
           gsl_vector_safe_free(u_hat);
           gsl_vector_safe_free(e_hat);
           gsl_vector_safe_free(y_hat);
-        }
+        } // output residuals
       }
 
       // Fit LMM or mvLMM (w. LOCO)
-      if (cPar.a_mode == 1 || cPar.a_mode == 2 || cPar.a_mode == 3 ||
-          cPar.a_mode == 4) {
+      if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 ||
+          cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM9) {
         if (cPar.n_ph == 1) {
+          debug_msg("fit LMM (one phenotype)");
           LMM cLmm;
-          cLmm.CopyFromParam(cPar);
+          cLmm.CopyFromParam(cPar); // set parameters
 
           // if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
@@ -2811,8 +2816,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
           cLmm.WriteFiles();
           cLmm.CopyToParam(cPar);
         } else {
+          debug_msg("fit mvLMM (multiple phenotypes)");
           MVLMM cMvlmm;
-          cMvlmm.CopyFromParam(cPar);
+          cMvlmm.CopyFromParam(cPar); // set parameters
 
           // if (is_check_mode()) disable_segfpe(); // disable fast NaN checking
 
diff --git a/src/gemma.h b/src/gemma.h
index 67a7f58..b34a958 100644
--- a/src/gemma.h
+++ b/src/gemma.h
@@ -41,8 +41,10 @@ using namespace std;
 // gw:      72
 
 enum M_MODE { M_LMM1=1, M_LMM2=2, M_LMM3=3, M_LMM4=4, M_LMM5=5,
-              M_BSLMM5=15,
-              M_KIN=21, M_KIN2=22, M_EIGEN=31 };
+  M_LMM9=9,  // GeneNetwork mode
+  M_BSLMM5=15,
+  M_KIN=21, M_KIN2=22, M_EIGEN=31
+};
 
 class GEMMA {
 
diff --git a/src/lmm.cpp b/src/lmm.cpp
index e10c0e9..9dd4aa6 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -42,6 +42,7 @@
 #include "gsl/gsl_vector.h"
 
 #include "gzstream.h"
+#include "gemma.h"
 #include "gemma_io.h"
 #include "fastblas.h"
 #include "lapack.h"
@@ -109,27 +110,27 @@ void LMM::WriteFiles() {
   }
 
   auto common_header = [&] () {
-    if (a_mode != 2) {
+    if (a_mode != M_LMM2) {
       outfile << "beta" << "\t";
       outfile << "se" << "\t";
     }
 
-    if (a_mode != 3)
-      outfile << "logl_H1" << "\t";  // we may make this an option
+    if (a_mode != M_LMM3)
+      outfile << "logl_H1" << "\t";
 
     switch(a_mode) {
-    case 1:
+    case M_LMM1:
       outfile << "l_remle" << "\t"
               << "p_wald" << endl;
       break;
-    case 2:
+    case M_LMM2:
       outfile << "l_mle" << "\t"
               << "p_lrt" << endl;
       break;
-    case 3:
+    case M_LMM3:
       outfile << "p_score" << endl;
       break;
-    case 4:
+    case M_LL4:
       outfile << "l_remle" << "\t"
               << "l_mle" << "\t"
               << "p_wald" << "\t"
@@ -1848,17 +1849,17 @@ void LMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
         FUNC_PARAM param1 = {false, ni_test, n_cvt, eval, Uab, ab, 0};
 
         // 3 is before 1, for beta.
-        if (a_mode == 3 || a_mode == 4) {
+        if (a_mode == M_LMM3 || a_mode == M_LMM4) {
           CalcRLScore(l_mle_null, param1, beta, se, p_score);
         }
 
-        if (a_mode == 1 || a_mode == 4) {
+        if (a_mode == M_LMM1 || a_mode == M_LMM4) {
           CalcLambda('R', param1, l_min, l_max, n_region, lambda_remle,
                      logl_H1);
           CalcRLWald(lambda_remle, param1, beta, se, p_wald);
         }
 
-        if (a_mode == 2 || a_mode == 4) {
+        if (a_mode == M_LMM2 || a_mode == M_LMM4) {
           CalcLambda('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
           p_lrt = gsl_cdf_chisq_Q(2.0 * (logl_H1 - logl_mle_H0), 1);
         }