gemma,lmm,mvlmm: Add data debug printing everywhere.

3 files changed, 19 insertions, 4 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index 4dfd7fd..48e377d 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1777,6 +1777,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     // write(eval,"eval zeroed");
     cPar.time_eigen = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
 
+    write(W, "W before UtW");
     // calculate UtW and Uty
     CalcUtX(U, W, UtW);
     CalcUtX(U, Y, UtY);
@@ -2552,7 +2553,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_LMM9 ||
-      cPar.a_mode == M_EIGEN  ) { // Fit LMM or mvLMM or eigen
+      cPar.a_mode == M_EIGEN) { // Fit LMM or mvLMM or eigen
+    write(cPar.a_mode, "Mode");
     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);
@@ -2589,6 +2591,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       // center matrix G
       CenterMatrix(G);
       validate_K(G);
+      write(G, "G");
 
       // is residual weights are provided, then
       if (!cPar.file_weight.empty()) {
@@ -2685,8 +2688,11 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cLmm.WriteFiles();
       cLmm.CopyToParam(cPar);
     } else {
-      debug_msg("Main LMM track");
+      write("Main LMM track", "Main LMM track");
       // calculate UtW and Uty
+      write(W, "W");
+      write(Y, "Y");
+      write(U, "U");
       CalcUtX(U, W, UtW);
       CalcUtX(U, Y, UtY);
       assert_issue(is_issue(26), ROUND(UtY->data[0]) == -16.6143);
@@ -2779,6 +2785,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
         } // output residuals
       }
 
+      write(UtW, "UtW");
+
       // Fit LMM or mvLMM (w. LOCO)
       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) {
@@ -2791,6 +2799,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
           gsl_vector_view Y_col = gsl_matrix_column(Y, 0);
           gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
+          write(&Y_col.vector, "Y_col");
 
           if (!cPar.file_bfile.empty()) {
             // PLINK analysis
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 091b3b4..85e92fe 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1478,6 +1478,8 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
                   const set<string> gwasnps) {
   clock_t time_start = clock();
 
+  write(W, "W");
+  write(y, "y");
   // Subset/LOCO support
   bool process_gwasnps = gwasnps.size();
   if (process_gwasnps)
@@ -2162,7 +2164,7 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
 
   gsl_matrix_set_zero(Uab);
   write(UtW,"UtW");
-  write(UtW,"Uty");
+  write(Uty,"Uty");
   CalcUab(UtW, Uty, Uab);
   write(Uab,"Uab");
   Calcab(UtW, Uty, ab);
@@ -2211,6 +2213,7 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
   size_t n_cvt = UtW->size2, ni_test = UtW->size1;
   size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
 
+  write(Uty, "VgVe Uty");
   gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index);
   gsl_vector *ab = gsl_vector_safe_alloc(n_index);
   gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index);
@@ -2238,7 +2241,8 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW,
   }
   fast_dgemm("T", "N", 1.0, HiW, UtW, 0.0, WHiW);
   gsl_blas_dgemv(CblasTrans, 1.0, HiW, Uty, 0.0, WHiy);
-
+  write(WHiW, "VgVe WHiW");
+  write(WHiy, "VgVe WHiy");
   int sig;
   gsl_permutation *pmt = gsl_permutation_alloc(UtW->size2);
   LUDecomp(WHiW, pmt, &sig);
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index 1d05bf1..f9cb0a2 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -2768,6 +2768,7 @@ void MphInitial(const size_t em_iter, const double em_prec,
                 gsl_matrix *B) {
 
   debug_msg("MphInitial");
+  write(Y, "Y in MphInitial");
   gsl_matrix_set_zero(V_g);
   gsl_matrix_set_zero(V_e);
   gsl_matrix_set_zero(B);
@@ -2971,6 +2972,7 @@ double PCRT(const size_t mode, const size_t d_size, const double p_value,
 void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval,
                           const gsl_matrix *UtW, const gsl_matrix *UtY) {
   debug_msg("entering");
+  write(UtY, "UtY in AnalyzeBimbam");
   igzstream infile(file_geno.c_str(), igzstream::in);
   if (!infile) {
     cout << "error reading genotype file:" << file_geno << endl;