diff options
| -rw-r--r-- | src/gemma_io.cpp | 2 | ||||
| -rw-r--r-- | src/lmm.cpp | 13 | ||||
| -rw-r--r-- | src/mathfunc.cpp | 6 |
3 files changed, 15 insertions, 6 deletions
diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp index 569c79b..81182f8 100644 --- a/src/gemma_io.cpp +++ b/src/gemma_io.cpp @@ -1557,6 +1557,7 @@ bool BimbamKin(const string file_geno, const set<string> ksnps, // FIXME: the following is not so slow but appears to generate an // identical matrix + /* for (size_t i = 0; i < ni_total; ++i) { for (size_t j = 0; j < i; ++j) { double d = gsl_matrix_get(matrix_kin, j, i); @@ -1564,6 +1565,7 @@ bool BimbamKin(const string file_geno, const set<string> ksnps, } } write(matrix_kin,"K rotated"); + */ // GSL is faster - and there are even faster methods // enforce_gsl(gsl_matrix_transpose(matrix_kin)); diff --git a/src/lmm.cpp b/src/lmm.cpp index 6337116..5e53fa2 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -277,9 +277,9 @@ 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, const gsl_matrix *Uab, const gsl_vector *unused, gsl_matrix *Pab) { - - // size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2; // result size - // auto ni_test = Uab->size1; // inds +#if !defined NDEBUG + 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); @@ -287,6 +287,7 @@ void CalcPab(const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval, assert(Pab->size2 == n_index); assert(Hi_eval->size == ni_test); // assert(ab->size == n_index); +#endif // DEBUG // compute Hi_eval (inds) * Uab (inds x n_index) * ab (n_index) and return in Pab (cvt x n_index). @@ -593,8 +594,9 @@ $7 = 3 $8 = 6 */ - // auto Uab = p->Uab; - // auto ab = p->ab; +#if !defined NDEBUG + 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? @@ -606,6 +608,7 @@ $8 = 6 assert(p->e_mode == 0); assert(Hi_eval->size == ni_test); +#endif // DEBUG 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); diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp index aaa9431..e74a841 100644 --- a/src/mathfunc.cpp +++ b/src/mathfunc.cpp @@ -400,7 +400,11 @@ uint count_abs_small_values(const gsl_vector *v, double min) { // and the ratio of max and min but one (min is expected to be zero). bool isMatrixIllConditioned(const gsl_vector *eigenvalues, double max_ratio) { auto t = abs_minmax(eigenvalues); - // auto absmin = get<0>(t); + +#if !defined NDEBUG + auto absmin = get<0>(t); +#endif + auto absmin1 = get<1>(t); auto absmax = get<2>(t); if (absmax/absmin1 > max_ratio) { |