diff options
| -rw-r--r-- | src/lmm.cpp | 16 |
1 files changed, 9 insertions, 7 deletions
diff --git a/src/lmm.cpp b/src/lmm.cpp index f849ee1..3f222d1 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -300,21 +300,23 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *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 + for (size_t p = 0; p <= n_cvt + 1; ++p) { // p walks rows phenotypes + covariates + for (size_t a = p + 1; a <= n_cvt + 2; ++a) { // a walks cols in p+1..rest for (size_t b = a; b <= n_cvt + 2; ++b) { // b in a..rest 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) { + if (p == 0) { // fills row 0 for each a,b using dot product of Hi_eval . Uab(a) + cout << "p is 0 " << index_ab; // walk row 0 + gsl_vector_const_view Uab_col = gsl_matrix_const_column(Uab, index_ab); // get the column + gsl_blas_ddot(Hi_eval, &Uab_col.vector, &p_ab); // dot product with H_eval + if (e_mode != 0) { // if not null model (defunct right now) if (! is_legacy_mode()) assert(false); // disabling to see when it is used; allow with legacy mode p_ab = gsl_vector_get(unused, index_ab) - p_ab; // was ab } cout << p << "r," << index_ab << ":" << p_ab << endl; gsl_matrix_set(Pab, 0, index_ab, p_ab); + write(Pab,"Pab int"); } else { + // walk the rest of the upper triangle of the matrix (row 1..n). Cols jump with 2 at a time cout << "a" << a << "b" << b << "p" << p << "n_cvt" << n_cvt << endl; write(Pab,"Pab int"); size_t index_aw = GetabIndex(a, p, n_cvt); |