Debugging calcPab

5 files changed, 194 insertions, 40 deletions

diff --git a/src/debug.cpp b/src/debug.cpp
index 37e4e8a..8df4334 100644
--- a/src/debug.cpp
+++ b/src/debug.cpp
@@ -148,6 +148,31 @@ void disable_segfpe() {
   #endif
 }
 
+void write(const gsl_vector *v, const char *msg) {
+  if (msg) cout << "// " << msg << endl;
+  cout << "// vector size: " << v->size << endl;
+  cout << "double " << msg << "[] = {";
+  for (size_t i=0; i < v->size; i++) {
+    cout << gsl_vector_get(v,i) << ",";
+  }
+  cout << "}" << endl;
+}
+
+void write(const gsl_matrix *m, const char *msg) {
+  if (msg) cout << "// " << msg << endl;
+  cout << "// matrix size: " << m->size1 << " cols, " << m->size2 << " rows" << endl;
+  cout << "double " << msg << "[] = {";
+  for (size_t j=0; j < m->size2; j++) {
+    for (size_t i=0; i < m->size1; i++) {
+      // cout << "(" << i << "," << j << ")";
+      cout << gsl_matrix_safe_get(m,i,j);
+      cout << ",";
+    }
+    cout << "// row " << j << endl;
+  }
+  cout << "}" << endl;
+}
+
 /*
   Helper function to make sure gsl allocations do their job because
   gsl_matrix_alloc does not initiatize values (behaviour that changed
diff --git a/src/debug.h b/src/debug.h
index 64ab721..a4be9f5 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -52,6 +52,9 @@ void disable_segfpe();
   { auto x = m * n;                                      \
     enforce_msg(x / m == n, "multiply integer overflow"); }
 
+void write(const gsl_vector *v, const char *msg);
+void write(const gsl_matrix *m, const char *msg);
+
 gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols);
 int gsl_matrix_safe_memcpy (gsl_matrix *dest, const gsl_matrix *src);
 void gsl_matrix_safe_free (gsl_matrix *v);
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 6c04c8d..ddb937e 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -544,9 +544,9 @@ void GEMMA::PrintHelp(size_t option) {
     cout << " -lmax     [num]          "
          << " specify maximum value for lambda (default 1e+5)" << endl;
     cout
-        << " -region    [num]          "
-        << " specify the number of regions used to evaluate lambda (default 10)"
-        << endl;
+         << " -region   [num]          "
+         << " specify the number of regions used to evaluate lambda (default 10)"
+         << endl;
     cout << " -loco     [chr]          "
          << " leave one chromosome out (LOCO) by name (requires -a annotation "
             "file)"
diff --git a/src/lmm.cpp b/src/lmm.cpp
index cdddf40..bba2e28 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -218,42 +218,119 @@ void LMM::WriteFiles() {
   return;
 }
 
+/*
+   As explained in
+   https://github.com/genetics-statistics/GEMMA/issues/94 CalcPab
+   returns the Pab matrix. As described For pab, it stores all
+   variables in the form of v_a P_p v_b. (Similarly, ppab stores all
+   v_a P_p P_p v_b, while pppab stores all v_a P_p P_p P_p v_b. These
+   quantities are defined according to the page 6 of this
+   supplementary information
+   http://xzlab.org/papers/2012_Zhou&Stephens_NG_SI.pdf).
+
+   In the code, p, a, b are indexes: when p=n_cvt+1, P_p is P_x as in
+   that supplementary information; when a=n_cvt+1, v_a=x; and when
+   a=n_cvt+2, v_a=y.
+
+   e_mode determines which model the algorithm is fitting: when
+   e_mode==1, it computes all the above quantities for the alternative
+   model (with the random effects term); when e_mode==0, it compute
+   these quantities for the null model (without the random effects
+   term).  Note that e==0 is only used here.
+
+   These quantities were computed based on the initial GEMMA paper,
+   and the goal is to finally compute y P_x y, y P_x P_xy,
+   y P_x P_x P_x y and the few trace forms (section 3.1.4 on page 5 of
+   the supplementary information). Sometimes I was wondering if we
+   should compute all these final quantities directly, instead of
+   through these complicated recursions. Direct computation may only
+   make computation a little slower, but will make the code much
+   easier to follow and easier to modify
+
+   a typical call sends n_cvt a vector
+   Hi_eval, a vector ab and a matrix Uab.
+
+  CalcPab(n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
+
+(gdb) p Uab->size1
+$1 = 247
+(gdb) p n_cvt
+$2 = 1
+(gdb) p e_mode
+$3 = 0
+(gdb) p Uab->size2
+$4 = 6
+(gdb) p Hi_eval->size
+$5 = 247
+(gdb) p ab->size
+$6 = 6
+(gdb) p Pab->size1
+$7 = 3
+(gdb) p Pab->size2
+$8 = 6
+
+  Hi_eval [0..ind] x Uab [ind, n_index] x ab [n_index]
+ */
+
 void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval,
              const gsl_matrix *Uab, const gsl_vector *ab, gsl_matrix *Pab) {
-  size_t index_ab, index_aw, index_bw, index_ww;
-  double p_ab;
-  double ps_ab, ps_aw, ps_bw, ps_ww;
+
+
+  size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2; // result size
+  auto ni_test = Uab->size1; // inds
+  assert(Uab->size1 == Hi_eval->size);
+  assert(Uab->size2 == n_index);
+
+  assert(Pab->size1 == n_cvt+2);
+  assert(Pab->size2 == n_index);
+  assert(Hi_eval->size == ni_test);
+  assert(ab->size == n_index);
+
+  // compute Hi_eval (inds)  * Uab  (inds x n_index) * ab (n_index) and return in Pab (cvt x n_index).
+
+  double p_ab = 0.0;
+
+  write(Hi_eval,"Hi_eval");
+  write(Uab,"Uab");
+  write(ab,"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
       for (size_t b = a; b <= n_cvt + 2; ++b) {    // b in a..rest
-        index_ab = GetabIndex(a, b, n_cvt);
+        size_t index_ab = GetabIndex(a, b, n_cvt); // index in top half matrix, see above
         if (p == 0) {
           gsl_vector_const_view Uab_col =
               gsl_matrix_const_column(Uab, index_ab);
           gsl_blas_ddot(Hi_eval, &Uab_col.vector, &p_ab);
           if (e_mode != 0) {
-            assert(false);
+            if (! is_legacy_mode()) assert(false); // disabling to see when it is used
             p_ab = gsl_vector_get(ab, index_ab) - p_ab;
           }
           gsl_matrix_set(Pab, 0, index_ab, p_ab);
         } else {
-          index_aw = GetabIndex(a, p, n_cvt);
-          index_bw = GetabIndex(b, p, n_cvt);
-          index_ww = GetabIndex(p, p, n_cvt);
-
-          ps_ab = gsl_matrix_safe_get(Pab, p - 1, index_ab);
-          ps_aw = gsl_matrix_safe_get(Pab, p - 1, index_aw);
-          ps_bw = gsl_matrix_safe_get(Pab, p - 1, index_bw);
-          ps_ww = gsl_matrix_safe_get(Pab, p - 1, index_ww);
+          size_t index_aw = GetabIndex(a, p, n_cvt);
+          size_t index_bw = GetabIndex(b, p, n_cvt);
+          size_t index_ww = GetabIndex(p, p, n_cvt);
+
+          auto rows = Pab->size1; // n_cvt+2
+          double ps_ab = 0.0;
+          if (index_ab < rows) ps_ab = gsl_matrix_safe_get(Pab, p - 1, index_ab);
+          double ps_aw = 0.0;
+          if (index_aw < rows) ps_aw = gsl_matrix_safe_get(Pab, p - 1, index_aw);
+          double ps_bw = 0.0;
+          if (index_bw < rows) ps_bw = gsl_matrix_safe_get(Pab, p - 1, index_bw);
+          double ps_ww = 0.0;
+          if (index_ww < rows) ps_ww = gsl_matrix_safe_get(Pab, p - 1, index_ww);
+
+          if (ps_ww != 0.0)
+            p_ab = ps_ab - ps_aw * ps_bw / ps_ww;
 
-          assert(ps_ww != 0.0);
-          p_ab = ps_ab - ps_aw * ps_bw / ps_ww;
           gsl_matrix_set(Pab, p, index_ab, p_ab);
         }
       }
     }
   }
+  write(Pab,"Pab");
   return;
 }
 
@@ -290,7 +367,6 @@ void CalcPPab(const size_t n_cvt, const size_t e_mode,
           ps2_bw = gsl_matrix_get(PPab, p - 1, index_bw);
           ps2_ww = gsl_matrix_get(PPab, p - 1, index_ww);
 
-          assert(ps_ww != 0.0);
           p2_ab = ps2_ab + ps_aw * ps_bw * ps2_ww / (ps_ww * ps_ww);
           p2_ab -= (ps_aw * ps2_bw + ps_bw * ps2_aw) / ps_ww;
           gsl_matrix_set(PPab, p, index_ab, p2_ab);
@@ -340,7 +416,6 @@ void CalcPPPab(const size_t n_cvt, const size_t e_mode,
           ps3_bw = gsl_matrix_get(PPPab, p - 1, index_bw);
           ps3_ww = gsl_matrix_get(PPPab, p - 1, index_ww);
 
-          assert(ps_ww != 0.0);
           p3_ab = ps3_ab -
                   ps_aw * ps_bw * ps2_ww * ps2_ww / (ps_ww * ps_ww * ps_ww);
           p3_ab -= (ps_aw * ps3_bw + ps_bw * ps3_aw + ps2_aw * ps2_bw) / ps_ww;
@@ -399,9 +474,15 @@ double LogL_f(double l, void *params) {
   index_yy = GetabIndex(n_cvt + 2, n_cvt + 2, n_cvt);
   double P_yy = gsl_matrix_get(Pab, nc_total, index_yy);
 
-  assert(!is_nan(P_yy));
+  if (is_check_mode() || is_debug_mode()) {
+    cerr << "P_yy is" << P_yy << endl;
+    assert(!is_nan(P_yy));
+    assert(P_yy > 0);
+  }
   f = c - 0.5 * logdet_h - 0.5 * (double)ni_test * log(P_yy);
-  assert(!is_nan(f));
+  if (is_check_mode() || is_debug_mode()) {
+    assert(!is_nan(f));
+  }
   gsl_matrix_safe_free(Pab); // FIXME
   gsl_vector_safe_free(Hi_eval);
   gsl_vector_safe_free(v_temp);
@@ -412,7 +493,7 @@ double LogL_dev1(double l, void *params) {
   FUNC_PARAM *p = (FUNC_PARAM *)params;
   size_t n_cvt = p->n_cvt;
   size_t ni_test = p->ni_test;
-  size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2;
+  size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2; // represents top half of covariate matrix
 
   size_t nc_total;
   if (p->calc_null == true) {
@@ -450,6 +531,39 @@ double LogL_dev1(double l, void *params) {
     trace_Hi = (double)ni_test - trace_Hi;
   }
 
+  /*
+(gdb) p Uab->size1
+$1 = 247
+(gdb) p n_cvt
+$2 = 1
+(gdb) p e_mode
+$3 = 0
+(gdb) p Uab->size2
+$4 = 6
+(gdb) p Hi_eval->size
+$5 = 247
+(gdb) p ab->size
+$6 = 6
+(gdb) p Pab->size1
+$7 = 3
+(gdb) p Pab->size2
+$8 = 6
+  */
+
+  auto Uab = p->Uab;
+  auto ab = p->ab;
+  assert(n_index == (n_cvt + 2 + 1) * (n_cvt + 2) / 2);
+  assert(Uab->size1 == ni_test);
+  assert(Uab->size2 == n_index); // n_cvt == 1 -> n_index == 6?
+
+  assert(Pab->size1 == n_cvt+2);
+  assert(Pab->size2 == n_index);
+
+  assert(ab->size == n_index);
+
+  assert(p->e_mode == 0);
+  assert(Hi_eval->size == ni_test);
+
   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);
 
@@ -1102,6 +1216,7 @@ 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;
@@ -1173,6 +1288,7 @@ 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,
diff --git a/src/param.cpp b/src/param.cpp
index 1a2b57c..6cbeab7 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -1353,27 +1353,37 @@ void compAKtoS(const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt,
   return;
 }
 
+/*
 // Copied from lmm.cpp; is used in the following function compKtoV
-// map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1
+// map a number 1..(n_cvt+2) to an index between 0 and [(n_cvt+2)*2+(n_cvt+2)]/2-1
+// or 1..cols to 0..(cols*2+cols)/2-1.
+
+  For a 3x3 matrix the following index gets returned to CalcPab:
+
+  CalcPab
+  * 1,1:0
+  * 1,2:1
+  * 1,3:2
+  * 2,2:3
+  * 2,3:4
+  * 3,3:5
+
+which is really the iteration moving forward along the diagonal and
+items to the right of it.
+  */
+
+
 size_t GetabIndex(const size_t a, const size_t b, const size_t n_cvt) {
-  if (a > n_cvt + 2 || b > n_cvt + 2 || a <= 0 || b <= 0) {
-    cout << "error in GetabIndex." << endl;
-    return 0;
-  }
-  size_t index;
-  size_t l, h;
-  if (b > a) {
-    l = a;
-    h = b;
-  } else {
-    l = b;
-    h = a;
+  auto cols = n_cvt + 2;
+  enforce_msg(a<=cols && b<=cols,"GetabIndex problem");
+  size_t a1 = a, b1 = b;
+  if (b <= a) {
+    a1 = b;
+    b1 = a;
   }
 
-  size_t n = n_cvt + 2;
-
-  assert(2+h-1 != 0);
-  index = (2 * n - l + 2) * (l - 1) / 2 + h - l;
+  size_t index = (2 * cols - a1 + 2) * (a1 - 1) / 2 + b1 - a1;
+  // cerr << "* " << a1 << "," << b1 << "," << cols << ":" << index << endl;
   return index;
 
   // return ( b < a ?  ((2 * n - b + 2) * (b - 1) / 2 + a - b ): ((2 * n - a + 2) * (a - 1) / 2 + b - a) );