diff options
| -rw-r--r-- | src/lapack.cpp | 9 | ||||
| -rw-r--r-- | src/mvlmm.cpp | 10 | ||||
| -rwxr-xr-x | test/dev_test_suite.sh | 8 |
3 files changed, 20 insertions, 7 deletions
diff --git a/src/lapack.cpp b/src/lapack.cpp index eb5b16b..31b6add 100644 --- a/src/lapack.cpp +++ b/src/lapack.cpp @@ -253,8 +253,10 @@ double EigenDecomp(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, return d; } -// Same as EigenDecomp but zeroes eigenvalues close to zero. When -// negative eigenvalues remain a warning is issued. +// Does NOT set eigenvalues to be positive. G gets destroyed. Returns +// eigen trace and values in U and eval (eigenvalues). Same as +// EigenDecomp but zeroes eigenvalues close to zero. When negative +// eigenvalues remain a warning is issued. double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t flag_largematrix) { EigenDecomp(G,U,eval,flag_largematrix); @@ -262,7 +264,8 @@ double EigenDecomp_Zeroed(gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, int count_zero_eigenvalues = 0; int count_negative_eigenvalues = 0; for (size_t i = 0; i < eval->size; i++) { - if (std::abs(gsl_vector_get(eval, i)) < EIGEN_MINVALUE) + // if (std::abs(gsl_vector_get(eval, i)) < EIGEN_MINVALUE) + if (gsl_vector_get(eval, i) < 1e-10) gsl_vector_set(eval, i, 0.0); // checks if (gsl_vector_get(eval,i) == 0.0) diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp index a3c7c85..4b8db05 100644 --- a/src/mvlmm.cpp +++ b/src/mvlmm.cpp @@ -256,7 +256,15 @@ double EigenProc(const gsl_matrix *V_g, const gsl_matrix *V_e, gsl_vector *D_l, gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, V_e_hi, VgVehi, 0.0, Lambda); // Eigen decomposition of Lambda. - EigenDecomp_Zeroed(Lambda, U_l, D_l, 0); + EigenDecomp(Lambda, U_l, D_l, 0); + + for (size_t i = 0; i < d_size; i++) { + d = gsl_vector_get(D_l, i); + if (d < 0) { + gsl_vector_set(D_l, i, 0); + } + } + // Calculate UltVeh and UltVehi. gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, U_l, V_e_h, 0.0, UltVeh); diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh index 0992635..34a1576 100755 --- a/test/dev_test_suite.sh +++ b/test/dev_test_suite.sh @@ -158,15 +158,17 @@ testPlinkUnivariateLinearMixedModelLOCO1() { testCorrelatedPhenotypesMvLLM() { + # https://github.com/genetics-statistics/GEMMA/issues/179 + outn=corrpheno $gemma $gemmaopts -p data/correlated_phenotypes/Ysim_reg_gemma.txt \ -g data/correlated_phenotypes/Genotypes_gemma.csv \ -d data/correlated_phenotypes/Kinship_eigenval_gemma.txt \ -u data/correlated_phenotypes/Kinship_eigenvec_gemma.txt \ - -lmm 2 -n 1 9 4 6 10 -o corrpheno + -lmm 2 -n 1 9 4 6 10 -o $outn assertEquals 0 $? - # outfn=output/$outn.assoc.txt + outfn=output/$outn.assoc.txt # assertEquals "68" `wc -l < $outfn` - # assertEquals "15465346.22" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn` + assertEquals "777.32" `perl -nle 'foreach $x (split(/\s+/,$_)) { $sum += sprintf("%.2f",(substr($x,,0,6))) } END { printf "%.2f",$sum }' $outfn` } |