diff options
| -rw-r--r-- | src/debug.cpp | 32 | ||||
| -rw-r--r-- | src/debug.h | 11 | ||||
| -rw-r--r-- | src/fastblas.cpp | 64 | ||||
| -rw-r--r-- | src/gemma.cpp | 36 | ||||
| -rw-r--r-- | src/io.cpp | 27 | ||||
| -rw-r--r-- | src/io.h | 6 | ||||
| -rw-r--r-- | src/lmm.h | 2 | ||||
| -rw-r--r-- | src/main.cpp | 2 | ||||
| -rw-r--r-- | src/mathfunc.cpp | 9 | ||||
| -rw-r--r-- | src/mvlmm.cpp | 7 | ||||
| -rw-r--r-- | src/param.cpp | 22 | ||||
| -rw-r--r-- | src/param.h | 14 | ||||
| -rwxr-xr-x | test/dev_test_suite.sh | 6 | ||||
| -rw-r--r-- | test/src/unittests-math.cpp | 1 |
14 files changed, 132 insertions, 107 deletions
diff --git a/src/debug.cpp b/src/debug.cpp index dacb89d..eb5c041 100644 --- a/src/debug.cpp +++ b/src/debug.cpp @@ -19,36 +19,46 @@ #include "mathfunc.h" static bool debug_mode = false; -static bool debug_no_check = false; -static bool debug_strict = false; +static bool debug_check = false; // check data/algorithms +static bool debug_strict = false; // fail on error +static bool debug_quiet = false; +static uint debug_issue = 0; // github issues +static bool debug_legacy = false; // legacy mode void debug_set_debug_mode(bool setting) { debug_mode = setting; } -void debug_set_no_check_mode(bool setting) {debug_no_check = setting; } +void debug_set_no_check_mode(bool setting) {debug_check = !setting; } void debug_set_strict_mode(bool setting) { debug_strict = setting; } +void debug_set_quiet_mode(bool setting) { debug_quiet = setting; } +void debug_set_issue(uint issue) { debug_issue = issue; } +void debug_set_legacy_mode(bool setting) { debug_legacy = setting; } bool is_debug_mode() { return debug_mode; }; -bool is_no_check_mode() { return debug_no_check; }; +bool is_no_check_mode() { return !debug_check; }; +bool is_check_mode() { return debug_check; }; bool is_strict_mode() { return debug_strict; }; +bool is_quiet_mode() { return debug_quiet; }; +bool is_issue(uint issue) { return issue == debug_issue; }; +bool is_legacy_mode() { return debug_legacy; }; /* Helper function to make sure gsl allocations do their job because gsl_matrix_alloc does not initiatize values (behaviour that changed in GSL2) we introduced a 'strict mode' by initializing the buffer - with NaNs. This happens in STRICT mode without NO-CHECKS - (i.e. -strict option). + with NaNs. This happens when NO-CHECKS is not set + (i.e. -no-checks option). */ gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols) { gsl_matrix *m = gsl_matrix_alloc(rows,cols); enforce_msg(m,"Not enough memory"); // just to be sure when there is no error handler set - if (debug_strict && !debug_no_check) { + if (is_check_mode()) { gsl_matrix_set_all(m, nan("")); } return m; } // Helper function called by macro validate_K(K, check) -void do_validate_K(const gsl_matrix *K, bool do_check, bool strict, const char *__file, int __line) { - if (do_check) { +void do_validate_K(const gsl_matrix *K, const char *__file, int __line) { + if (is_check_mode()) { // debug_msg("Validating K"); auto eigenvalues = getEigenValues(K); const uint count_small = count_small_values(eigenvalues,EIGEN_MINVALUE); @@ -61,13 +71,13 @@ void do_validate_K(const gsl_matrix *K, bool do_check, bool strict, const char * if (!isMatrixIllConditioned(eigenvalues)) warning_at_msg(__file,__line,"K is ill conditioned!"); if (!isMatrixSymmetric(K)) - fail_at_msg(strict,__file,__line,"K is not symmetric!" ); + fail_at_msg(is_strict_mode(),__file,__line,"K is not symmetric!" ); const bool negative_values = has_negative_values_but_one(eigenvalues); if (negative_values) { warning_at_msg(__file,__line,"K has more than one negative eigenvalues!"); } if (count_small>1 && negative_values && !isMatrixPositiveDefinite(K)) - fail_at_msg(strict,__file,__line,"K is not positive definite!"); + fail_at_msg(is_strict_mode(),__file,__line,"K is not positive definite!"); gsl_vector_free(eigenvalues); } } diff --git a/src/debug.h b/src/debug.h index 0f7b38f..ee36887 100644 --- a/src/debug.h +++ b/src/debug.h @@ -13,18 +13,25 @@ void gemma_gsl_error_handler (const char * reason, void debug_set_debug_mode(bool setting); void debug_set_no_check_mode(bool setting); void debug_set_strict_mode(bool setting); +void debug_set_quiet_mode(bool setting); +void debug_set_issue(uint issue); +void debug_set_legacy_mode(bool setting); bool is_debug_mode(); bool is_no_check_mode(); +bool is_check_mode(); bool is_strict_mode(); +bool is_quiet_mode(); +bool is_issue(uint issue); +bool is_legacy_mode(); gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols); // Validation routines -void do_validate_K(const gsl_matrix *K, bool do_check, bool strict, const char *__file, int __line); +void do_validate_K(const gsl_matrix *K, const char *__file, int __line); #define ROUND(f) round(f * 10000.)/10000 -#define validate_K(K,check,strict) do_validate_K(K,check,strict,__FILE__,__LINE__) +#define validate_K(K) do_validate_K(K,__FILE__,__LINE__) #define warning_at_msg(__file,__line,msg) cerr << "**** WARNING: " << msg << " in " << __file << " at line " << __line << endl; diff --git a/src/fastblas.cpp b/src/fastblas.cpp index b25d964..0a4eba3 100644 --- a/src/fastblas.cpp +++ b/src/fastblas.cpp @@ -44,14 +44,16 @@ gsl_matrix *fast_copy(gsl_matrix *m, const double *mem) { } } } else { // faster goes by row + auto v = gsl_vector_calloc(cols); + enforce(v); // just to be sure for (auto r=0; r<rows; r++) { - auto v = gsl_vector_calloc(cols); assert(v->size == cols); assert(v->block->size == cols); assert(v->stride == 1); memcpy(v->block->data,&mem[r*cols],cols*sizeof(double)); gsl_matrix_set_row(m,r,v); } + gsl_vector_free(v); } return m; } @@ -160,28 +162,25 @@ static void fast_cblas_dgemm(const char *TransA, const char *TransB, const doubl // C++ is row-major auto transA = (*TransA == 'N' || *TransA == 'n' ? CblasNoTrans : CblasTrans); auto transB = (*TransB == 'N' || *TransB == 'n' ? CblasNoTrans : CblasTrans); - // A(m x k) * B(k x n) = C(m x n)) - auto rowsA = A->size1; - auto colsA = A->size2; - blasint M = A->size1; - blasint K = B->size1; - assert(K == colsA); - blasint N = B->size2; - // cout << M << "," << N "," << K << endl; - // Layout = CblasRowMajor: Trans: K , NoTrans M - blasint lda = (transA==CblasNoTrans ? K : M ); - blasint ldb = (transB==CblasNoTrans ? N : K ); - blasint ldc = N; - - fast_cblas_dgemm(CblasRowMajor, transA, transB, M, N, K, alpha, - /* A */ A->data, - /* lda */ lda, - /* B */ B->data, - /* ldb */ ldb, - /* beta */ beta, - /* C */ C->data, ldc); + const size_t M = C->size1; + const size_t N = C->size2; + const size_t MA = (transA == CblasNoTrans) ? A->size1 : A->size2; + const size_t NA = (transA == CblasNoTrans) ? A->size2 : A->size1; + const size_t MB = (transB == CblasNoTrans) ? B->size1 : B->size2; + const size_t NB = (transB == CblasNoTrans) ? B->size2 : B->size1; + + if (M == MA && N == NB && NA == MB) { /* [MxN] = [MAxNA][MBxNB] */ + + cblas_dgemm (CblasRowMajor, transA, transB, M, N,NA, + alpha, A->data, A->tda, B->data, B->tda, beta, + C->data, C->tda); + + } else { + throw invalid_argument("Range error in dgemm"); + } } + /* Use the fasted/supported way to call BLAS dgemm */ @@ -189,14 +188,19 @@ static void fast_cblas_dgemm(const char *TransA, const char *TransB, const doubl void fast_dgemm(const char *TransA, const char *TransB, const double alpha, const gsl_matrix *A, const gsl_matrix *B, const double beta, gsl_matrix *C) { - fast_cblas_dgemm(TransA,TransB,alpha,A,B,beta,C); - - #ifndef NDEBUG - if (is_strict_mode() && !is_no_check_mode()) { - // ---- validate with original implementation - gsl_matrix *C1 = gsl_matrix_alloc(C->size1,C->size2); - eigenlib_dgemm(TransA,TransB,alpha,A,B,beta,C1); - enforce_msg(gsl_matrix_equal(C,C1),"dgemm outcomes are not equal for fast & eigenlib"); + if (is_legacy_mode()) { + eigenlib_dgemm(TransA,TransB,alpha,A,B,beta,C); + } else { + fast_cblas_dgemm(TransA,TransB,alpha,A,B,beta,C); + + #ifndef NDEBUG + if (is_check_mode()) { + // ---- validate with original implementation + gsl_matrix *C1 = gsl_matrix_alloc(C->size1,C->size2); + eigenlib_dgemm(TransA,TransB,alpha,A,B,beta,C1); + enforce_msg(gsl_matrix_equal(C,C1),"dgemm outcomes are not equal for fast & eigenlib"); + gsl_matrix_free(C1); + } + #endif } - #endif } diff --git a/src/gemma.cpp b/src/gemma.cpp index 95630c6..e2881a4 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -721,6 +721,7 @@ void GEMMA::PrintHelp(size_t option) { cout << " -debug debug output" << endl; cout << " -nind [num] read up to num individuals" << endl; cout << " -issue [num] enable tests relevant to issue tracker" << endl; + cout << " -legacy run gemma in legacy mode" << endl; cout << endl; } @@ -766,7 +767,7 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) { str.assign(argv[i]); cPar.file_mbfile = str; } else if (strcmp(argv[i], "-silence") == 0) { - cPar.mode_silence = true; + debug_set_quiet_mode(true); } else if (strcmp(argv[i], "-g") == 0) { if (argv[i + 1] == NULL || argv[i + 1][0] == '-') { continue; @@ -1367,8 +1368,9 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) { ++i; str.clear(); str.assign(argv[i]); - cPar.issue = atoi(str.c_str()); // for testing purposes - enforce(cPar.issue > 0); + auto issue = atoi(str.c_str()); // for testing purposes + enforce(issue > 0); + debug_set_issue(issue); } else if (strcmp(argv[i], "-emp") == 0) { if (argv[i + 1] == NULL || argv[i + 1][0] == '-') { continue; @@ -1588,14 +1590,16 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) { str.assign(argv[i]); cPar.window_ns = atoi(str.c_str()); } else if (strcmp(argv[i], "-debug") == 0) { - cPar.mode_debug = true; + // cPar.mode_debug = true; debug_set_debug_mode(true); } else if (strcmp(argv[i], "-no-check") == 0) { - cPar.mode_check = false; + // cPar.mode_check = false; debug_set_no_check_mode(true); } else if (strcmp(argv[i], "-strict") == 0) { - cPar.mode_strict = true; + // cPar.mode_strict = true; debug_set_strict_mode(true); + } else if (strcmp(argv[i], "-legacy") == 0) { + debug_set_legacy_mode(true); } else { cout << "error! unrecognized option: " << argv[i] << endl; cPar.error = true; @@ -1742,7 +1746,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // center matrix G CenterMatrix(G); CenterMatrix(G_full); - validate_K(G,cPar.mode_check,cPar.mode_strict); + validate_K(G); // eigen-decomposition and calculate trace_G cout << "Start Eigen-Decomposition..." << endl; @@ -1882,7 +1886,7 @@ void GEMMA::BatchRun(PARAM &cPar) { } // Now we have the Kinship matrix test it - validate_K(G,cPar.mode_check,cPar.mode_strict); + validate_K(G); if (cPar.a_mode == 21) { cPar.WriteMatrix(G, "cXX"); @@ -2323,7 +2327,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // center matrix G CenterMatrix(G); - validate_K(G,cPar.mode_check,cPar.mode_strict); + validate_K(G); (cPar.v_traceG).clear(); double d = 0; @@ -2532,7 +2536,7 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_vector *eval = gsl_vector_calloc(Y->size1); gsl_vector *env = gsl_vector_alloc(Y->size1); gsl_vector *weight = gsl_vector_alloc(Y->size1); - assert_issue(cPar.issue == 26, UtY->data[0] == 0.0); + assert_issue(is_issue(26), UtY->data[0] == 0.0); // set covariates matrix W and phenotype matrix Y // an intercept should be included in W, @@ -2552,7 +2556,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // center matrix G CenterMatrix(G); - validate_K(G,cPar.mode_check,cPar.mode_strict); + validate_K(G); // is residual weights are provided, then if (!cPar.file_weight.empty()) { @@ -2633,7 +2637,7 @@ void GEMMA::BatchRun(PARAM &cPar) { CalcUtX(U, W, UtW); CalcUtX(U, Y, UtY); - assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143); + assert_issue(is_issue(26), ROUND(UtY->data[0]) == -16.6143); LMM cLmm; cLmm.CopyFromParam(cPar); @@ -2650,7 +2654,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // calculate UtW and Uty CalcUtX(U, W, UtW); CalcUtX(U, Y, UtY); - assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143); + assert_issue(is_issue(26), ROUND(UtY->data[0]) == -16.6143); // calculate REMLE/MLE estimate and pve for univariate model if (cPar.n_ph == 1) { // one phenotype @@ -2658,7 +2662,7 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_vector_view se_beta = gsl_matrix_row(se_B, 0); gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); - assert_issue(cPar.issue == 26, ROUND(UtY->data[0]) == -16.6143); + assert_issue(is_issue(26), ROUND(UtY->data[0]) == -16.6143); CalcLambda('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0); @@ -2857,7 +2861,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // center matrix G CenterMatrix(G); - validate_K(G,cPar.mode_check,cPar.mode_strict); + validate_K(G); } else { cPar.ReadGenotypes(UtX, G, true); } @@ -2968,7 +2972,7 @@ void GEMMA::BatchRun(PARAM &cPar) { // center matrix G CenterMatrix(G); - validate_K(G,cPar.mode_check,cPar.mode_strict); + validate_K(G); } else { cPar.ReadGenotypes(UtX, G, true); @@ -41,6 +41,7 @@ #include "debug.h" #include "eigenlib.h" +#include "fastblas.h" #include "gzstream.h" #include "io.h" #include "lapack.h" @@ -594,7 +595,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps, const double &r2_level, map<string, string> &mapRS2chr, map<string, long int> &mapRS2bp, map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo, - size_t &ns_test, bool debug) { + size_t &ns_test) { debug_msg("entered"); indicator_snp.clear(); snpInfo.clear(); @@ -667,7 +668,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps, } if (mapRS2bp.count(rs) == 0) { - if (debug && count_warnings++ < 10) { + if (is_debug_mode() && count_warnings++ < 10) { std::string msg = "Can't figure out position for "; msg += rs; debug_msg(msg); @@ -1454,12 +1455,12 @@ bool BimbamKin(const string file_geno, const set<string> ksnps, // compute kinship matrix and return in matrix_kin a SNP at a time if (ns_test % msize == 0) { - eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); + fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); gsl_matrix_set_zero(Xlarge); } } if (ns_test % msize != 0) { - eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); + fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); } cout << endl; @@ -1600,13 +1601,13 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp, ns_test++; if (ns_test % msize == 0) { - eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); + fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); gsl_matrix_set_zero(Xlarge); } } if (ns_test % msize != 0) { - eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); + fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); } cout << endl; @@ -1633,7 +1634,7 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp, // genotype and calculate K. bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv, vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, - const bool calc_K, bool debug) { + const bool calc_K) { debug_msg("entered"); igzstream infile(file_geno.c_str(), igzstream::in); if (!infile) { @@ -1738,7 +1739,7 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<vector<unsigned char>> &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, - const size_t ns_test, bool debug) { + const size_t ns_test) { debug_msg("entered"); igzstream infile(file_geno.c_str(), igzstream::in); if (!infile) { @@ -2757,7 +2758,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps, ns_vec[0]++; if (ns_vec[0] % msize == 0) { - eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); + fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); gsl_matrix_set_zero(Xlarge); } } else if (mapRS2cat.count(rs) != 0) { @@ -2774,7 +2775,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps, gsl_matrix_submatrix(Xlarge, 0, msize * i_vc, ni_test, msize); gsl_matrix_view kin_sub = gsl_matrix_submatrix( matrix_kin, 0, ni_test * i_vc, ni_test, ni_test); - eigenlib_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, + fast_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, &kin_sub.matrix); gsl_matrix_set_zero(&X_sub.matrix); @@ -2790,7 +2791,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps, gsl_matrix_submatrix(Xlarge, 0, msize * i_vc, ni_test, msize); gsl_matrix_view kin_sub = gsl_matrix_submatrix(matrix_kin, 0, ni_test * i_vc, ni_test, ni_test); - eigenlib_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, + fast_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, &kin_sub.matrix); } } @@ -2983,7 +2984,7 @@ bool PlinkKin(const string &file_bed, const int display_pace, ns_vec[0]++; if (ns_vec[0] % msize == 0) { - eigenlib_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); + fast_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin); gsl_matrix_set_zero(Xlarge); } } else if (mapRS2cat.count(rs) != 0) { @@ -3000,7 +3001,7 @@ bool PlinkKin(const string &file_bed, const int display_pace, gsl_matrix_submatrix(Xlarge, 0, msize * i_vc, ni_test, msize); gsl_matrix_view kin_sub = gsl_matrix_submatrix( matrix_kin, 0, ni_test * i_vc, ni_test, ni_test); - eigenlib_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, + fast_dgemm("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, &kin_sub.matrix); gsl_matrix_set_zero(&X_sub.matrix); @@ -64,7 +64,7 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps, const double &r2_level, map<string, string> &mapRS2chr, map<string, long int> &mapRS2bp, map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo, - size_t &ns_test, bool debug); + size_t &ns_test); bool ReadFile_bed(const string &file_bed, const set<string> &setSnps, const gsl_matrix *W, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<SNPINFO> &snpInfo, @@ -94,7 +94,7 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp, bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv, vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, - const bool calc_K, bool debug); + const bool calc_K); bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv, vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, const bool calc_K); @@ -102,7 +102,7 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<vector<unsigned char>> &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, - const size_t ns_test, bool debug); + const size_t ns_test); bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<vector<unsigned char>> &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, @@ -28,7 +28,7 @@ using namespace std; -#define LMM_BATCH_SIZE 1000 // used for batch processing +#define LMM_BATCH_SIZE 20000 // used for batch processing class FUNC_PARAM { diff --git a/src/main.cpp b/src/main.cpp index 92c4d90..d752a72 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -61,7 +61,7 @@ int main(int argc, char *argv[]) { return EXIT_FAILURE; } - if (cPar.mode_silence) { + if (is_quiet_mode()) { stringstream ss; cout.rdbuf(ss.rdbuf()); } diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp index fbbc061..d0a9d65 100644 --- a/src/mathfunc.cpp +++ b/src/mathfunc.cpp @@ -49,6 +49,7 @@ #include "debug.h" #include "eigenlib.h" +#include "fastblas.h" #include "lapack.h" #include "mathfunc.h" @@ -388,20 +389,16 @@ void StandardizeVector(gsl_vector *y) { void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) { gsl_matrix *X = gsl_matrix_alloc(UtX->size1, UtX->size2); gsl_matrix_memcpy(X, UtX); - eigenlib_dgemm("T", "N", 1.0, U, X, 0.0, UtX); + fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX); gsl_matrix_free(X); - - return; } void CalcUtX(const gsl_matrix *U, const gsl_matrix *X, gsl_matrix *UtX) { - eigenlib_dgemm("T", "N", 1.0, U, X, 0.0, UtX); - return; + fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX); } void CalcUtX(const gsl_matrix *U, const gsl_vector *x, gsl_vector *Utx) { gsl_blas_dgemv(CblasTrans, 1.0, U, x, 0.0, Utx); - return; } // Kronecker product. diff --git a/src/mvlmm.cpp b/src/mvlmm.cpp index c0dc143..95ef14b 100644 --- a/src/mvlmm.cpp +++ b/src/mvlmm.cpp @@ -39,6 +39,7 @@ #include "gsl/gsl_vector.h" #include "eigenlib.h" +#include "fastblas.h" #include "gzstream.h" #include "io.h" #include "lapack.h" @@ -3250,8 +3251,8 @@ void MVLMM::AnalyzeBimbam(const gsl_matrix *U, const gsl_vector *eval, gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l); time_start = clock(); - eigenlib_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, - &UtXlarge_sub.matrix); + fast_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, + &UtXlarge_sub.matrix); time_UtX += (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0); gsl_matrix_set_zero(Xlarge); @@ -3717,7 +3718,7 @@ void MVLMM::AnalyzePlink(const gsl_matrix *U, const gsl_vector *eval, gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l); time_start = clock(); - eigenlib_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, + fast_dgemm("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, &UtXlarge_sub.matrix); time_UtX += (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0); diff --git a/src/param.cpp b/src/param.cpp index d637481..8aa7a64 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -66,7 +66,7 @@ void LOCO_set_Snps(set<string> &ksnps, set<string> &gwasnps, // (indicator_idv[x] == 1). This should match indicator_cvt etc. If // this gives problems with certain sets we can simply trim to size. -void trim_individuals(vector<int> &idvs, size_t ni_max, bool debug) { +void trim_individuals(vector<int> &idvs, size_t ni_max) { if (ni_max) { size_t count = 0; for (auto ind = idvs.begin(); ind != idvs.end(); ++ind) { @@ -76,7 +76,7 @@ void trim_individuals(vector<int> &idvs, size_t ni_max, bool debug) { break; } if (count != idvs.size()) { - if (debug) + if (is_debug_mode()) cout << "**** TEST MODE: trim individuals from " << idvs.size() << " to " << count << endl; idvs.resize(count); @@ -87,7 +87,7 @@ void trim_individuals(vector<int> &idvs, size_t ni_max, bool debug) { // ---- PARAM class implementation PARAM::PARAM(void) - : mode_silence(false), a_mode(0), k_mode(1), d_pace(DEFAULT_PACE), + : a_mode(0), k_mode(1), d_pace(DEFAULT_PACE), file_out("result"), path_out("./output/"), miss_level(0.05), maf_level(0.01), hwe_level(0), r2_level(0.9999), l_min(1e-5), l_max(1e5), n_region(10), p_nr(0.001), em_prec(0.0001), nr_prec(0.0001), @@ -221,7 +221,7 @@ void PARAM::ReadFiles(void) { } else { n_cvt = 1; } - trim_individuals(indicator_cvt, ni_max, mode_debug); + trim_individuals(indicator_cvt, ni_max); if (!file_gxe.empty()) { if (ReadFile_column(file_gxe, indicator_gxe, gxe, 1) == false) { @@ -234,7 +234,7 @@ void PARAM::ReadFiles(void) { } } - trim_individuals(indicator_idv, ni_max, mode_debug); + trim_individuals(indicator_idv, ni_max); // Read genotype and phenotype file for PLINK format. if (!file_bfile.empty()) { @@ -302,11 +302,11 @@ void PARAM::ReadFiles(void) { gsl_matrix *W = gsl_matrix_alloc(ni_test, n_cvt); CopyCvt(W); - trim_individuals(indicator_idv, ni_max, mode_debug); - trim_individuals(indicator_cvt, ni_max, mode_debug); + trim_individuals(indicator_idv, ni_max); + trim_individuals(indicator_cvt, ni_max); if (ReadFile_geno(file_geno, setSnps, W, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, - mapRS2bp, mapRS2cM, snpInfo, ns_test, mode_debug) == false) { + mapRS2bp, mapRS2cM, snpInfo, ns_test) == false) { error = true; } gsl_matrix_free(W); @@ -416,7 +416,7 @@ void PARAM::ReadFiles(void) { while (!safeGetline(infile, file_name).eof()) { if (ReadFile_geno(file_name, setSnps, W, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, - mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp, mode_debug) == false) { + mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp) == false) { error = true; } @@ -1282,7 +1282,7 @@ void PARAM::ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) { } } else { if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, UtX, K, - calc_K, mode_debug) == false) { + calc_K) == false) { error = true; } } @@ -1302,7 +1302,7 @@ void PARAM::ReadGenotypes(vector<vector<unsigned char>> &Xt, gsl_matrix *K, } } else { if (ReadFile_geno(file_geno, indicator_idv, indicator_snp, Xt, K, calc_K, - ni_test, ns_test, mode_debug) == false) { + ni_test, ns_test) == false) { error = true; } } diff --git a/src/param.h b/src/param.h index 4b473c0..c4316bb 100644 --- a/src/param.h +++ b/src/param.h @@ -26,8 +26,8 @@ #include <set> #include <vector> -#define K_BATCH_SIZE 10000 // #snps used for batched K -#define DEFAULT_PACE 1000 +#define K_BATCH_SIZE 20000 // #snps used for batched K +#define DEFAULT_PACE 1000 // for display only using namespace std; @@ -116,11 +116,11 @@ public: class PARAM { public: // IO-related parameters - bool mode_check = true; // run data checks (slower) - bool mode_strict = false; // exit on some data checks - bool mode_silence; - bool mode_debug = false; - uint issue; // enable tests for issue on github tracker + // bool mode_check = true; // run data checks (slower) + // bool mode_strict = false; // exit on some data checks + // bool mode_silence; + // bool mode_debug = false; + // uint issue; // enable tests for issue on github tracker uint a_mode; // Analysis mode, 1/2/3/4 for Frequentist tests int k_mode; // Kinship read mode: 1: n by n matrix, 2: id/id/k_value; diff --git a/test/dev_test_suite.sh b/test/dev_test_suite.sh index 0ccc24c..77604be 100755 --- a/test/dev_test_suite.sh +++ b/test/dev_test_suite.sh @@ -11,7 +11,7 @@ testBXDStandardRelatednessMatrixKSingularError() { -c ../example/BXD_covariates.txt \ -a ../example/BXD_snps.txt \ -gk \ - -debug -strict \ + -debug \ -o $outn assertEquals 22 $? # should show singular error } @@ -24,7 +24,7 @@ testBXDStandardRelatednessMatrixK() { -c ../example/BXD_covariates2.txt \ -a ../example/BXD_snps.txt \ -gk \ - -debug -strict \ + -debug \ -o $outn assertEquals 0 $? outfn=output/$outn.cXX.txt @@ -40,7 +40,7 @@ testBXDLMMLikelihoodRatio() { -a ../example/BXD_snps.txt \ -k ./output/BXD.cXX.txt \ -lmm 2 -maf 0.1 \ - -debug -strict \ + -debug \ -o $outn assertEquals 0 $? diff --git a/test/src/unittests-math.cpp b/test/src/unittests-math.cpp index 6486707..757c2dc 100644 --- a/test/src/unittests-math.cpp +++ b/test/src/unittests-math.cpp @@ -10,6 +10,7 @@ #include "debug.h" #include "mathfunc.h" #include "fastblas.h" +#include "fastopenblas.h" using namespace std; |