From 59f453afe05b0b637edeff3127fb7ef7be5944b4 Mon Sep 17 00:00:00 2001 From: Pjotr Prins Date: Fri, 31 Aug 2018 12:29:45 +0000 Subject: Improving readability --- src/gemma.cpp | 43 ++++++++++++++++++++++++++++++------------- src/lmm.cpp | 2 +- src/mathfunc.cpp | 2 +- src/mvlmm.cpp | 2 -- 4 files changed, 32 insertions(+), 17 deletions(-) diff --git a/src/gemma.cpp b/src/gemma.cpp index 94d05dc..e1d76d4 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -64,6 +64,23 @@ extern "C" { using namespace std; +// OPTIONS +// ------- +// gk: 21-22 +// gs: 25-26 +// gq: 27-28 +// eigen: 31-32 +// lmm: 1-5 +// bslmm: 11-15 +// predict: 41-43 +// lm: 51 +// vc: 61 +// ci: 66-67 +// calccor: 71 +// gw: 72 + +enum M_MODE { M_LMM1=1, M_LMM2=2, M_LMM3=3, M_LMM4=4, M_LMM5=5, M_KIN=21, M_KIN2=22, M_EIGEN=31 }; + GEMMA::GEMMA(void) : version(GEMMA_VERSION), date(GEMMA_DATE), year(GEMMA_YEAR) {} void gemma_gsl_error_handler (const char * reason, @@ -1126,7 +1143,7 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) { break; } if (argv[i + 1] == NULL || argv[i + 1][0] == '-') { - cPar.a_mode = 21; + cPar.a_mode = M_KIN; continue; } ++i; @@ -1878,7 +1895,7 @@ void GEMMA::BatchRun(PARAM &cPar) { } // Generate Kinship matrix (optionally using LOCO) - if (cPar.a_mode == 21 || cPar.a_mode == 22) { + if (cPar.a_mode == M_KIN || cPar.a_mode == M_KIN2) { cout << "Calculating Relatedness Matrix ... " << endl; gsl_matrix *G = gsl_matrix_safe_alloc(cPar.ni_total, cPar.ni_total); @@ -1897,7 +1914,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // Now we have the Kinship matrix test it validate_K(G); - if (cPar.a_mode == 21) { + if (cPar.a_mode == M_KIN) { cPar.WriteMatrix(G, "cXX"); } else { cPar.WriteMatrix(G, "sXX"); @@ -2535,9 +2552,9 @@ void GEMMA::BatchRun(PARAM &cPar) { } // LMM or mvLMM or Eigen-Decomposition - if (cPar.a_mode == 1 || cPar.a_mode == 2 || cPar.a_mode == 3 || - cPar.a_mode == 4 || cPar.a_mode == 5 || - cPar.a_mode == 31) { // Fit LMM or mvLMM or eigen + if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 || + cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 || + cPar.a_mode == M_EIGEN) { // Fit LMM or mvLMM or eigen gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph); enforce_msg(Y, "allocate Y"); // just to be sure gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt); @@ -2599,7 +2616,7 @@ void GEMMA::BatchRun(PARAM &cPar) { cout << "Start Eigen-Decomposition..." << endl; time_start = clock(); - if (cPar.a_mode == 31) { + if (cPar.a_mode == M_EIGEN) { cPar.trace_G = EigenDecomp_Zeroed(G, U, eval, 1); } else { cPar.trace_G = EigenDecomp_Zeroed(G, U, eval, 0); @@ -2646,7 +2663,7 @@ void GEMMA::BatchRun(PARAM &cPar) { } // write(eval,"eval2"); - if (cPar.a_mode == 31) { + if (cPar.a_mode == M_EIGEN) { cPar.WriteMatrix(U, "eigenU"); cPar.WriteVector(eval, "eigenD"); } else if (!cPar.file_gene.empty()) { @@ -2725,7 +2742,7 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.PrintSummary(); // calculate and output residuals - if (cPar.a_mode == 5) { + if (cPar.a_mode == M_LMM5) { gsl_vector *Utu_hat = gsl_vector_safe_alloc(Y->size1); gsl_vector *Ute_hat = gsl_vector_safe_alloc(Y->size1); gsl_vector *u_hat = gsl_vector_safe_alloc(Y->size1); @@ -3514,17 +3531,17 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) { outfile << "## total computation time = " << cPar.time_total << " min " << endl; outfile << "## computation time break down: " << endl; - if (cPar.a_mode == 21 || cPar.a_mode == 22 || cPar.a_mode == 11 || + if (cPar.a_mode == M_KIN || cPar.a_mode == M_KIN2 || cPar.a_mode == 11 || cPar.a_mode == 13 || cPar.a_mode == 14 || cPar.a_mode == 16) { outfile << "## time on calculating relatedness matrix = " << cPar.time_G << " min " << endl; } - if (cPar.a_mode == 31) { + if (cPar.a_mode == M_EIGEN) { outfile << "## time on eigen-decomposition = " << cPar.time_eigen << " min " << endl; } - if (cPar.a_mode == 1 || cPar.a_mode == 2 || cPar.a_mode == 3 || - cPar.a_mode == 4 || cPar.a_mode == 5 || cPar.a_mode == 11 || + if (cPar.a_mode == M_LMM1 || cPar.a_mode == M_LMM2 || cPar.a_mode == M_LMM3 || + cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 || cPar.a_mode == 11 || cPar.a_mode == 12 || cPar.a_mode == 13 || cPar.a_mode == 14 || cPar.a_mode == 16) { outfile << "## time on eigen-decomposition = " << cPar.time_eigen diff --git a/src/lmm.cpp b/src/lmm.cpp index b368c0b..6e14346 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -1955,7 +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); + // Calcab(W, y, ab); FUNC_PARAM param0 = {true, ni_test, n_cvt, eval, Uab, ab, 0}; diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp index 7117700..4708ffa 100644 --- a/src/mathfunc.cpp +++ b/src/mathfunc.cpp @@ -449,7 +449,7 @@ void StandardizeVector(gsl_vector *y) { return; } -// Calculate UtX. +// Calculate UtX (U gets transposed) void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) { gsl_matrix *X = gsl_matrix_safe_alloc(UtX->size1, UtX->size2); gsl_matrix_safe_memcpy(X, UtX); diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp index 0a259a7..13325ab 100644 --- a/src/mvlmm.cpp +++ b/src/mvlmm.cpp @@ -2817,8 +2817,6 @@ 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); -- cgit v1.2.3