Continue debugging calcpab

5 files changed, 73 insertions, 21 deletions

diff --git a/src/debug.cpp b/src/debug.cpp
index 8df4334..3efcce6 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -160,15 +160,22 @@ void write(const gsl_vector *v, const char *msg) {
 
 void write(const gsl_matrix *m, const char *msg) {
   if (msg) cout << "// " << msg << endl;
-  cout << "// matrix size: " << m->size1 << " cols, " << m->size2 << " rows" << endl;
+  // Matrices are stored in row-major order, meaning that each row of
+  // elements forms a contiguous block in memory. This is the standard
+  // “C-language ordering” of two-dimensional arrays. The number of
+  // rows is size1.
+  auto rows = m->size1; // see https://www.gnu.org/software/gsl/manual/html_node/Matrices.html#Matrices
+  auto cols = m->size2;
+
+  cout << "// matrix size: " << cols << " cols, " << rows << " rows" << endl;
   cout << "double " << msg << "[] = {";
-  for (size_t j=0; j < m->size2; j++) {
-    for (size_t i=0; i < m->size1; i++) {
+  for (size_t row=0; row < rows; row++) {
+    for (size_t col=0; col < cols; col++) {
       // cout << "(" << i << "," << j << ")";
-      cout << gsl_matrix_safe_get(m,i,j);
+      cout << gsl_matrix_safe_get(m,row,col);
       cout << ",";
     }
-    cout << "// row " << j << endl;
+    cout << "// row " << row << endl;
   }
   cout << "}" << endl;
 }
@@ -254,21 +261,21 @@ gsl_vector *gsl_vector_safe_alloc(size_t n) {
   return v;
 }
 
-double do_gsl_matrix_safe_get(const gsl_matrix * m, const size_t i, const size_t j,
+double do_gsl_matrix_safe_get(const gsl_matrix * m, const size_t row, const size_t col,
                               const char *__pretty_function, const char *__file, int __line) {
   enforce(m);
   if (!is_legacy_mode() && (is_debug_mode() || is_check_mode() || is_strict_mode())) {
-    auto rows = m->size1;
+    auto rows = m->size1; // see above write function
     auto cols = m->size2;
-    if (i >= cols || j >= rows) {
-      std::string msg = "Matrix out of bounds (" + std::to_string(cols) + "," + std::to_string(rows) + ") ";
-      msg += std::to_string(i);
-      msg += ",";
-      msg += std::to_string(j);
+    if (col >= cols || row >= rows) {
+      std::string msg = "Matrix out of bounds (" + std::to_string(rows) + "," + std::to_string(cols) + ") ";
+      msg += std::to_string(row);
+      msg += "r,";
+      msg += std::to_string(col);
       fail_at_msg(__file,__line,msg.c_str());
     }
   }
-  return gsl_matrix_get(m,i,j);
+  return gsl_matrix_get(m,row,col);
 }
 
 
diff --git a/src/gemma.cpp b/src/gemma.cpp
index ddb937e..94d05dc 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1759,6 +1759,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
     cout << "Start Eigen-Decomposition..." << endl;
     time_start = clock();
     cPar.trace_G = EigenDecomp_Zeroed(G, U, eval, 0);
+    write(eval,"eval zeroed");
     cPar.time_eigen = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
 
     // calculate UtW and Uty
@@ -2603,6 +2604,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       } else {
         cPar.trace_G = EigenDecomp_Zeroed(G, U, eval, 0);
       }
+      write(eval,"eval");
 
       if (!cPar.file_weight.empty()) {
         double wi;
@@ -2642,6 +2644,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
       }
       cPar.trace_G /= (double)eval->size;
     }
+    // write(eval,"eval2");
 
     if (cPar.a_mode == 31) {
       cPar.WriteMatrix(U, "eigenU");
@@ -2700,6 +2703,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
         assert(!std::isnan(cPar.beta_mle_null.front()));
         assert(!std::isnan(cPar.se_beta_mle_null.front()));
 
+        // the following functions do not modify eval
         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,
@@ -2798,6 +2802,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
           // if (is_check_mode()) disable_segfpe(); // disable fast NaN checking
 
+          // write(eval,"eval3");
+
           if (!cPar.file_bfile.empty()) {
             if (cPar.file_gxe.empty()) {
               cMvlmm.AnalyzePlink(U, eval, UtW, UtY);
diff --git a/src/lapack.cpp b/src/lapack.cpp
index 6e43bd9..165a82d 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -239,6 +239,8 @@ void lapack_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
 double EigenDecomp(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
                    const size_t flag_largematrix) {
   lapack_eigen_symmv(G, eval, U, flag_largematrix);
+  assert(!has_nan(eval));
+  // write(eval,"eval");
 
   // Calculate track_G=mean(diag(G)).
   double d = 0.0;
diff --git a/src/lmm.cpp b/src/lmm.cpp
index bba2e28..b368c0b 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -270,6 +270,8 @@ $7 = 3
 $8 = 6
 
   Hi_eval [0..ind] x Uab [ind, n_index] x ab [n_index]
+
+Iterating through a dataset Hi_eval differs and Uab (last row)
  */
 
 void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
@@ -290,9 +292,13 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
 
   double p_ab = 0.0;
 
+  debug_msg("CalcPab");
   write(Hi_eval,"Hi_eval");
   write(Uab,"Uab");
   write(ab,"ab");
+  assert(!has_nan(Hi_eval));
+  assert(!has_nan(Uab));
+  // assert(!has_nan(ab));
 
   for (size_t p = 0; p <= n_cvt + 1; ++p) {        // p walks phenotypes + covariates
     for (size_t a = p + 1; a <= n_cvt + 2; ++a) {  // a in p+1..rest
@@ -303,9 +309,10 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
               gsl_matrix_const_column(Uab, index_ab);
           gsl_blas_ddot(Hi_eval, &Uab_col.vector, &p_ab);
           if (e_mode != 0) {
-            if (! is_legacy_mode()) assert(false); // disabling to see when it is used
+            if (! is_legacy_mode()) assert(false); // disabling to see when it is used; allow with legacy mode
             p_ab = gsl_vector_get(ab, index_ab) - p_ab;
           }
+          cout << p << "r," << index_ab << ":" << p_ab << endl;
           gsl_matrix_set(Pab, 0, index_ab, p_ab);
         } else {
           size_t index_aw = GetabIndex(a, p, n_cvt);
@@ -325,12 +332,15 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
           if (ps_ww != 0.0)
             p_ab = ps_ab - ps_aw * ps_bw / ps_ww;
 
+          cout << p << "r," << index_ab << ":" << p_ab << endl;
           gsl_matrix_set(Pab, p, index_ab, p_ab);
         }
       }
     }
   }
   write(Pab,"Pab");
+  if (is_strict_mode() && (has_nan(Uab) || has_nan(Pab) || has_nan(Hi_eval) || has_nan(ab)))
+    exit(2);
   return;
 }
 
@@ -341,6 +351,11 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode,
   double p2_ab;
   double ps2_ab, ps_aw, ps_bw, ps_ww, ps2_aw, ps2_bw, ps2_ww;
 
+  debug_msg("CalcPPab");
+  write(HiHi_eval,"Hi_eval");
+  write(Uab,"Uab");
+  write(ab,"ab");
+
   for (size_t p = 0; p <= n_cvt + 1; ++p) {
     for (size_t a = p + 1; a <= n_cvt + 2; ++a) {
       for (size_t b = a; b <= n_cvt + 2; ++b) {
@@ -374,6 +389,10 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode,
       }
     }
   }
+  write(Pab,"Pab");
+  write(PPab,"PPab");
+  if (is_strict_mode() && (has_nan(Uab) || has_nan(PPab) || has_nan(HiHi_eval) || has_nan(ab)))
+    exit(2);
   return;
 }
 
@@ -385,6 +404,10 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode,
   double p3_ab;
   double ps3_ab, ps_aw, ps_bw, ps_ww, ps2_aw, ps2_bw, ps2_ww, ps3_aw, ps3_bw,
       ps3_ww;
+  debug_msg("CalcPPPab");
+  write(HiHiHi_eval,"HiHiHi_eval");
+  write(Uab,"Uab");
+  write(ab,"ab");
 
   for (size_t p = 0; p <= n_cvt + 1; ++p) {
     for (size_t a = p + 1; a <= n_cvt + 2; ++a) {
@@ -428,6 +451,11 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode,
       }
     }
   }
+  write(Pab,"Pab");
+  write(PPab,"PPab");
+  write(PPPab,"PPPab");
+  if (is_strict_mode() && (has_nan(Uab) || has_nan(PPPab) || has_nan(HiHiHi_eval) || has_nan(ab)))
+    exit(2);
   return;
 }
 
@@ -832,7 +860,8 @@ double LogRL_dev1(double l, void *params) {
   gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size);
   gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size);
 
-  gsl_vector_safe_memcpy(v_temp, p->eval);
+  // write(p->eval, "p->eval");
+  gsl_vector_safe_memcpy(v_temp, p->eval); // initialize with eval
   gsl_vector_scale(v_temp, l);
   if (p->e_mode == 0) {
     gsl_vector_set_all(Hi_eval, 1.0);
@@ -852,6 +881,8 @@ double LogRL_dev1(double l, void *params) {
     trace_Hi = (double)ni_test - trace_Hi;
   }
 
+  write(p->eval, "p->eval2");
+  write(p->ab, "p->ab");
   CalcPab(n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
   CalcPPab(n_cvt, p->e_mode, HiHi_eval, p->Uab, p->ab, Pab, PPab);
 
@@ -1216,12 +1247,9 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty,
   return;
 }
 
-/*
 void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
-  size_t index_ab;
   size_t n_cvt = W->size2;
 
-  double d;
   gsl_vector *v_a = gsl_vector_safe_alloc(y->size);
   gsl_vector *v_b = gsl_vector_safe_alloc(y->size);
 
@@ -1242,7 +1270,7 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
         continue;
       }
 
-      index_ab = GetabIndex(a, b, n_cvt);
+      auto index_ab = GetabIndex(a, b, n_cvt);
 
       if (b == n_cvt + 2) {
         gsl_vector_safe_memcpy(v_b, y);
@@ -1251,6 +1279,7 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
         gsl_vector_safe_memcpy(v_b, &W_col.vector);
       }
 
+      double d;
       gsl_blas_ddot(v_a, v_b, &d);
       gsl_vector_set(ab, index_ab, d);
     }
@@ -1288,7 +1317,6 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x,
   gsl_vector_safe_free(v_b);
   return;
 }
-*/
 
 void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval,
                       const gsl_matrix *UtW, const gsl_vector *Utx,
@@ -1766,7 +1794,7 @@ void CalcLambda(const char func_name, FUNC_PARAM &params, const double l_min,
     lambda_l = l_min * exp(lambda_interval * i);
     lambda_h = l_min * exp(lambda_interval * (i + 1.0));
 
-    if (func_name == 'R' || func_name == 'r') {
+    if (func_name == 'R' || func_name == 'r') { // log-restricted likelihood
       dev1_l = LogRL_dev1(lambda_l, &params);
       dev1_h = LogRL_dev1(lambda_h, &params);
     } else {
@@ -1909,6 +1937,9 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
                 const gsl_matrix *UtW, const gsl_vector *Uty,
                 const double l_min, const double l_max, const size_t n_region,
                 double &lambda, double &logl_H0) {
+
+  write(eval,"eval6");
+
   if (func_name != 'R' && func_name != 'L' && func_name != 'r' &&
       func_name != 'l') {
     cout << "func_name only takes 'R' or 'L': 'R' for "
@@ -1924,6 +1955,7 @@ void CalcLambda(const char func_name, const gsl_vector *eval,
 
   gsl_matrix_set_zero(Uab);
   CalcUab(UtW, Uty, Uab);
+  Calcab(UtW, Uty, ab);
 
   FUNC_PARAM param0 = {true, ni_test, n_cvt, eval, Uab, ab, 0};
 
diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp
index 14f8b4a..0a259a7 100644
--- a/src/mvlmm.cpp
+++ b/src/mvlmm.cpp
@@ -2754,6 +2754,7 @@ void MphInitial(const size_t em_iter, const double em_prec,
                 const size_t n_region, gsl_matrix *V_g, gsl_matrix *V_e,
                 gsl_matrix *B) {
 
+  debug_msg("MphInitial");
   gsl_matrix_set_zero(V_g);
   gsl_matrix_set_zero(V_e);
   gsl_matrix_set_zero(B);
@@ -2816,6 +2817,8 @@ void MphInitial(const size_t em_iter, const double em_prec,
     gsl_matrix *Ve_sub = gsl_matrix_alloc(2, 2);
     gsl_matrix *B_sub = gsl_matrix_alloc(2, c_size);
 
+    write(eval,"eval5");
+
     for (size_t i = 0; i < d_size; i++) {
       gsl_vector_view Y_sub1 = gsl_matrix_row(Y_sub, 0);
       gsl_vector_const_view Y_1 = gsl_matrix_const_row(Y, i);
@@ -3477,6 +3480,8 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval,
   MphInitial(em_iter, em_prec, nr_iter, nr_prec, eval, &X_sub.matrix, Y, l_min,
              l_max, n_region, V_g, V_e, &B_sub.matrix);
 
+  write(eval,"eval4");
+
   logl_H0 = MphEM('R', em_iter, em_prec, eval, &X_sub.matrix, Y, U_hat, E_hat,
                   OmegaU, OmegaE, UltVehiY, UltVehiBX, UltVehiU, UltVehiE, V_g,
                   V_e, &B_sub.matrix);