diff options
Diffstat (limited to 'src/gemma.cpp')
| -rw-r--r-- | src/gemma.cpp | 303 |
1 files changed, 158 insertions, 145 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp index 24173c3..f9a2fc9 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -23,6 +23,17 @@ #include <iostream> #include <string> #include <sys/stat.h> +#ifdef OPENBLAS +#pragma message "Compiling with OPENBLAS" +extern "C" { + // these functions are defined in cblas.h - but if we include that we + // conflicts with other BLAS includes + int openblas_get_num_threads(void); + int openblas_get_parallel(void); + char* openblas_get_config(void); + char* openblas_get_corename(void); +} +#endif #include "gsl/gsl_blas.h" #include "gsl/gsl_cdf.h" @@ -310,11 +321,6 @@ void GEMMA::PrintHelp(size_t option) { cout << " rs#2, base_position, chr_number" << endl; cout << " ..." << endl; - // WJA added. - cout << " -oxford [prefix] " - << " specify input Oxford genotype bgen file prefix." << endl; - cout << " requires: *.bgen, *.sample files" << endl; - cout << " -gxe [filename] " << " specify input file that contains a column of environmental " "factor for g by e tests" @@ -429,8 +435,8 @@ void GEMMA::PrintHelp(size_t option) { "default 1)" << endl; cout << " -pace [num] " - << " specify terminal display update pace (default 100000 SNPs or " - "100000 iterations)." + << " specify terminal display update pace (default 1,000 SNPs or " + "1,000 iterations)." << endl; cout << " -outdir [path] " << " specify output directory path (default \"./output/\")" << endl; @@ -715,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; } @@ -760,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; @@ -793,18 +800,6 @@ void GEMMA::Assign(int argc, char **argv, PARAM &cPar) { str.clear(); str.assign(argv[i]); cPar.file_anno = str; - } - - // WJA added. - else if (strcmp(argv[i], "-oxford") == 0 || - strcmp(argv[i], "--oxford") == 0 || strcmp(argv[i], "-x") == 0) { - if (argv[i + 1] == NULL || argv[i + 1][0] == '-') { - continue; - } - ++i; - str.clear(); - str.assign(argv[i]); - cPar.file_oxford = str; } else if (strcmp(argv[i], "-gxe") == 0) { if (argv[i + 1] == NULL || argv[i + 1][0] == '-') { continue; @@ -1373,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; @@ -1594,11 +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; @@ -1635,7 +1636,7 @@ void GEMMA::BatchRun(PARAM &cPar) { if (cPar.a_mode == 41 || cPar.a_mode == 42) { gsl_vector *y_prdt; - y_prdt = gsl_vector_alloc(cPar.ni_total - cPar.ni_test); + y_prdt = gsl_vector_safe_alloc(cPar.ni_total - cPar.ni_test); // set to zero gsl_vector_set_zero(y_prdt); @@ -1647,8 +1648,8 @@ void GEMMA::BatchRun(PARAM &cPar) { if (!cPar.file_kin.empty() && !cPar.file_ebv.empty()) { cout << "Adding Breeding Values ... " << endl; - gsl_matrix *G = gsl_matrix_alloc(cPar.ni_total, cPar.ni_total); - gsl_vector *u_hat = gsl_vector_alloc(cPar.ni_test); + gsl_matrix *G = gsl_matrix_safe_alloc(cPar.ni_total, cPar.ni_total); + gsl_vector *u_hat = gsl_vector_safe_alloc(cPar.ni_test); // read kinship matrix and set u_hat vector<int> indicator_all; @@ -1706,25 +1707,25 @@ void GEMMA::BatchRun(PARAM &cPar) { if (cPar.a_mode == 43) { // first, use individuals with full phenotypes to obtain estimates of Vg and // Ve - gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph); - gsl_matrix *W = gsl_matrix_alloc(Y->size1, cPar.n_cvt); - gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1); - gsl_matrix *U = gsl_matrix_alloc(Y->size1, Y->size1); - gsl_matrix *UtW = gsl_matrix_alloc(Y->size1, W->size2); - gsl_matrix *UtY = gsl_matrix_alloc(Y->size1, Y->size2); - gsl_vector *eval = gsl_vector_alloc(Y->size1); + gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph); + gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt); + gsl_matrix *G = gsl_matrix_safe_alloc(Y->size1, Y->size1); + gsl_matrix *U = gsl_matrix_safe_alloc(Y->size1, Y->size1); + gsl_matrix *UtW = gsl_matrix_safe_alloc(Y->size1, W->size2); + gsl_matrix *UtY = gsl_matrix_safe_alloc(Y->size1, Y->size2); + gsl_vector *eval = gsl_vector_safe_alloc(Y->size1); - gsl_matrix *Y_full = gsl_matrix_alloc(cPar.ni_cvt, cPar.n_ph); - gsl_matrix *W_full = gsl_matrix_alloc(Y_full->size1, cPar.n_cvt); + gsl_matrix *Y_full = gsl_matrix_safe_alloc(cPar.ni_cvt, cPar.n_ph); + gsl_matrix *W_full = gsl_matrix_safe_alloc(Y_full->size1, cPar.n_cvt); // set covariates matrix W and phenotype matrix Y // an intercept should be included in W, cPar.CopyCvtPhen(W, Y, 0); cPar.CopyCvtPhen(W_full, Y_full, 1); - gsl_matrix *Y_hat = gsl_matrix_alloc(Y_full->size1, cPar.n_ph); - gsl_matrix *G_full = gsl_matrix_alloc(Y_full->size1, Y_full->size1); - gsl_matrix *H_full = gsl_matrix_alloc(Y_full->size1 * Y_hat->size2, + gsl_matrix *Y_hat = gsl_matrix_safe_alloc(Y_full->size1, cPar.n_ph); + gsl_matrix *G_full = gsl_matrix_safe_alloc(Y_full->size1, Y_full->size1); + gsl_matrix *H_full = gsl_matrix_safe_alloc(Y_full->size1 * Y_hat->size2, Y_full->size1 * Y_hat->size2); // read relatedness matrix G, and matrix G_full @@ -1745,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; @@ -1760,8 +1761,8 @@ void GEMMA::BatchRun(PARAM &cPar) { // calculate variance component and beta estimates // and then obtain predicted values if (cPar.n_ph == 1) { - gsl_vector *beta = gsl_vector_alloc(W->size2); - gsl_vector *se_beta = gsl_vector_alloc(W->size2); + gsl_vector *beta = gsl_vector_safe_alloc(W->size2); + gsl_vector *se_beta = gsl_vector_safe_alloc(W->size2); double lambda, logl, vg, ve; gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); @@ -1791,10 +1792,10 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_vector_free(beta); gsl_vector_free(se_beta); } else { - gsl_matrix *Vg = gsl_matrix_alloc(cPar.n_ph, cPar.n_ph); - gsl_matrix *Ve = gsl_matrix_alloc(cPar.n_ph, cPar.n_ph); - gsl_matrix *B = gsl_matrix_alloc(cPar.n_ph, W->size2); - gsl_matrix *se_B = gsl_matrix_alloc(cPar.n_ph, W->size2); + gsl_matrix *Vg = gsl_matrix_safe_alloc(cPar.n_ph, cPar.n_ph); + gsl_matrix *Ve = gsl_matrix_safe_alloc(cPar.n_ph, cPar.n_ph); + gsl_matrix *B = gsl_matrix_safe_alloc(cPar.n_ph, W->size2); + gsl_matrix *se_B = gsl_matrix_safe_alloc(cPar.n_ph, W->size2); // obtain estimates CalcMvLmmVgVeBeta(eval, UtW, UtY, cPar.em_iter, cPar.nr_iter, @@ -1872,7 +1873,7 @@ void GEMMA::BatchRun(PARAM &cPar) { if (cPar.a_mode == 21 || cPar.a_mode == 22) { cout << "Calculating Relatedness Matrix ... " << endl; - gsl_matrix *G = gsl_matrix_alloc(cPar.ni_total, cPar.ni_total); + gsl_matrix *G = gsl_matrix_safe_alloc(cPar.ni_total, cPar.ni_total); enforce_msg(G, "allocate G"); // just to be sure time_start = clock(); @@ -1885,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"); @@ -1917,8 +1918,8 @@ void GEMMA::BatchRun(PARAM &cPar) { if (cPar.a_mode == 25 || cPar.a_mode == 26) { cout << "Calculating the S Matrix ... " << endl; - gsl_matrix *S = gsl_matrix_alloc(cPar.n_vc * 2, cPar.n_vc); - gsl_vector *ns = gsl_vector_alloc(cPar.n_vc + 1); + gsl_matrix *S = gsl_matrix_safe_alloc(cPar.n_vc * 2, cPar.n_vc); + gsl_vector *ns = gsl_vector_safe_alloc(cPar.n_vc + 1); gsl_matrix_set_zero(S); gsl_vector_set_zero(ns); @@ -1927,13 +1928,13 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_matrix_submatrix(S, cPar.n_vc, 0, cPar.n_vc, cPar.n_vc); gsl_vector_view ns_vec = gsl_vector_subvector(ns, 0, cPar.n_vc); - gsl_matrix *K = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); - gsl_matrix *A = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); + gsl_matrix *K = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); + gsl_matrix *A = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); gsl_matrix_set_zero(K); gsl_matrix_set_zero(A); - gsl_vector *y = gsl_vector_alloc(cPar.ni_test); - gsl_matrix *W = gsl_matrix_alloc(cPar.ni_test, cPar.n_cvt); + gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test); + gsl_matrix *W = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_cvt); cPar.CopyCvtPhen(W, y, 0); @@ -1970,9 +1971,9 @@ void GEMMA::BatchRun(PARAM &cPar) { // Compute the q vector, that is used for variance component estimation using // summary statistics if (cPar.a_mode == 27 || cPar.a_mode == 28) { - gsl_matrix *Vq = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc); - gsl_vector *q = gsl_vector_alloc(cPar.n_vc); - gsl_vector *s = gsl_vector_alloc(cPar.n_vc + 1); + gsl_matrix *Vq = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc); + gsl_vector *q = gsl_vector_safe_alloc(cPar.n_vc); + gsl_vector *s = gsl_vector_safe_alloc(cPar.n_vc + 1); gsl_vector_set_zero(q); gsl_vector_set_zero(s); @@ -2028,8 +2029,8 @@ void GEMMA::BatchRun(PARAM &cPar) { // LM. if (cPar.a_mode == 51 || cPar.a_mode == 52 || cPar.a_mode == 53 || cPar.a_mode == 54) { // Fit LM - gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph); - gsl_matrix *W = gsl_matrix_alloc(Y->size1, cPar.n_cvt); + gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph); + gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt); // set covariates matrix W and phenotype matrix Y // an intercept should be included in W, @@ -2047,8 +2048,6 @@ void GEMMA::BatchRun(PARAM &cPar) { &Y_col.vector); // y is the predictor, not the phenotype } else if (!cPar.file_bfile.empty()) { cLm.AnalyzePlink(W, &Y_col.vector); - } else if (!cPar.file_oxford.empty()) { - cLm.Analyzebgen(W, &Y_col.vector); } else { cLm.AnalyzeBimbam(W, &Y_col.vector); } @@ -2083,16 +2082,16 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.UpdateSNP(mapRS2wK); // Setup matrices and vectors. - gsl_matrix *S = gsl_matrix_alloc(cPar.n_vc * 2, cPar.n_vc); - gsl_matrix *Vq = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc); - gsl_vector *q = gsl_vector_alloc(cPar.n_vc); - gsl_vector *s = gsl_vector_alloc(cPar.n_vc + 1); + gsl_matrix *S = gsl_matrix_safe_alloc(cPar.n_vc * 2, cPar.n_vc); + gsl_matrix *Vq = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc); + gsl_vector *q = gsl_vector_safe_alloc(cPar.n_vc); + gsl_vector *s = gsl_vector_safe_alloc(cPar.n_vc + 1); - gsl_matrix *K = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); - gsl_matrix *A = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); + gsl_matrix *K = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); + gsl_matrix *A = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc * cPar.ni_test); - gsl_vector *y = gsl_vector_alloc(cPar.ni_test); - gsl_matrix *W = gsl_matrix_alloc(cPar.ni_test, cPar.n_cvt); + gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test); + gsl_matrix *W = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_cvt); gsl_matrix_set_zero(K); gsl_matrix_set_zero(A); @@ -2219,16 +2218,16 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.n_vc = cPar.n_vc - 1; - gsl_matrix *S = gsl_matrix_alloc(2 * cPar.n_vc, cPar.n_vc); - gsl_matrix *Vq = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc); - // gsl_matrix *V=gsl_matrix_alloc (cPar.n_vc+1, + gsl_matrix *S = gsl_matrix_safe_alloc(2 * cPar.n_vc, cPar.n_vc); + gsl_matrix *Vq = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc); + // gsl_matrix *V=gsl_matrix_safe_alloc (cPar.n_vc+1, // (cPar.n_vc*(cPar.n_vc+1))/2*(cPar.n_vc+1) ); - // gsl_matrix *Vslope=gsl_matrix_alloc (n_lines+1, + // gsl_matrix *Vslope=gsl_matrix_safe_alloc (n_lines+1, // (n_lines*(n_lines+1))/2*(n_lines+1) ); - gsl_vector *q = gsl_vector_alloc(cPar.n_vc); - gsl_vector *s_study = gsl_vector_alloc(cPar.n_vc); - gsl_vector *s_ref = gsl_vector_alloc(cPar.n_vc); - gsl_vector *s = gsl_vector_alloc(cPar.n_vc + 1); + gsl_vector *q = gsl_vector_safe_alloc(cPar.n_vc); + gsl_vector *s_study = gsl_vector_safe_alloc(cPar.n_vc); + gsl_vector *s_ref = gsl_vector_safe_alloc(cPar.n_vc); + gsl_vector *s = gsl_vector_safe_alloc(cPar.n_vc + 1); gsl_matrix_set_zero(S); gsl_matrix_view S_mat = @@ -2287,9 +2286,9 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_vector_free(s_ref); gsl_vector_free(s); } else { - gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph); - gsl_matrix *W = gsl_matrix_alloc(Y->size1, cPar.n_cvt); - gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1 * cPar.n_vc); + gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph); + gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt); + gsl_matrix *G = gsl_matrix_safe_alloc(Y->size1, Y->size1 * cPar.n_vc); // set covariates matrix W and phenotype matrix Y // an intercept should be included in W, @@ -2328,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; @@ -2366,9 +2365,9 @@ void GEMMA::BatchRun(PARAM &cPar) { // the genotypes if (cPar.a_mode == 66 || cPar.a_mode == 67) { // read reference file first - gsl_matrix *S = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc); - gsl_matrix *Svar = gsl_matrix_alloc(cPar.n_vc, cPar.n_vc); - gsl_vector *s_ref = gsl_vector_alloc(cPar.n_vc); + gsl_matrix *S = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc); + gsl_matrix *Svar = gsl_matrix_safe_alloc(cPar.n_vc, cPar.n_vc); + gsl_vector *s_ref = gsl_vector_safe_alloc(cPar.n_vc); gsl_matrix_set_zero(S); gsl_matrix_set_zero(Svar); @@ -2393,14 +2392,14 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.ObtainWeight(setSnps_beta, mapRS2wK); // set up matrices and vector - gsl_matrix *Xz = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc); - gsl_matrix *XWz = gsl_matrix_alloc(cPar.ni_test, cPar.n_vc); + gsl_matrix *Xz = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc); + gsl_matrix *XWz = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_vc); gsl_matrix *XtXWz = - gsl_matrix_alloc(mapRS2wK.size(), cPar.n_vc * cPar.n_vc); - gsl_vector *w = gsl_vector_alloc(mapRS2wK.size()); - gsl_vector *w1 = gsl_vector_alloc(mapRS2wK.size()); - gsl_vector *z = gsl_vector_alloc(mapRS2wK.size()); - gsl_vector *s_vec = gsl_vector_alloc(cPar.n_vc); + gsl_matrix_safe_alloc(mapRS2wK.size(), cPar.n_vc * cPar.n_vc); + gsl_vector *w = gsl_vector_safe_alloc(mapRS2wK.size()); + gsl_vector *w1 = gsl_vector_safe_alloc(mapRS2wK.size()); + gsl_vector *z = gsl_vector_safe_alloc(mapRS2wK.size()); + gsl_vector *s_vec = gsl_vector_safe_alloc(cPar.n_vc); vector<size_t> vec_cat, vec_size; vector<double> vec_z; @@ -2524,20 +2523,20 @@ void GEMMA::BatchRun(PARAM &cPar) { 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 - gsl_matrix *Y = gsl_matrix_alloc(cPar.ni_test, cPar.n_ph); + 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_alloc(Y->size1, cPar.n_cvt); - gsl_matrix *B = gsl_matrix_alloc(Y->size2, W->size2); // B is a d by c + gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt); + gsl_matrix *B = gsl_matrix_safe_alloc(Y->size2, W->size2); // B is a d by c // matrix - gsl_matrix *se_B = gsl_matrix_alloc(Y->size2, W->size2); - gsl_matrix *G = gsl_matrix_alloc(Y->size1, Y->size1); - gsl_matrix *U = gsl_matrix_alloc(Y->size1, Y->size1); + gsl_matrix *se_B = gsl_matrix_safe_alloc(Y->size2, W->size2); + gsl_matrix *G = gsl_matrix_safe_alloc(Y->size1, Y->size1); + gsl_matrix *U = gsl_matrix_safe_alloc(Y->size1, Y->size1); gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2); gsl_matrix *UtY = gsl_matrix_calloc(Y->size1, Y->size2); 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); + gsl_vector *env = gsl_vector_safe_alloc(Y->size1); + gsl_vector *weight = gsl_vector_safe_alloc(Y->size1); + 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, @@ -2557,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()) { @@ -2638,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); @@ -2655,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 @@ -2663,31 +2662,27 @@ 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); assert(!std::isnan(UtY->data[0])); - assert(!std::isnan(B->data[0])); - assert(!std::isnan(se_B->data[0])); CalcLmmVgVeBeta(eval, UtW, &UtY_col.vector, cPar.l_mle_null, cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector, &se_beta.vector); assert(!std::isnan(UtY->data[0])); - assert(!std::isnan(B->data[0])); - assert(!std::isnan(se_B->data[0])); cPar.beta_mle_null.clear(); cPar.se_beta_mle_null.clear(); + assert(!std::isnan(B->data[0])); + assert(!std::isnan(se_B->data[0])); for (size_t i = 0; i < B->size2; i++) { cPar.beta_mle_null.push_back(gsl_matrix_get(B, 0, i)); cPar.se_beta_mle_null.push_back(gsl_matrix_get(se_B, 0, i)); } assert(!std::isnan(UtY->data[0])); - assert(!std::isnan(B->data[0])); - assert(!std::isnan(se_B->data[0])); assert(!std::isnan(cPar.beta_mle_null.front())); assert(!std::isnan(cPar.se_beta_mle_null.front())); @@ -2699,6 +2694,9 @@ void GEMMA::BatchRun(PARAM &cPar) { cPar.beta_remle_null.clear(); cPar.se_beta_remle_null.clear(); + assert(!std::isnan(B->data[0])); + assert(!std::isnan(se_B->data[0])); + for (size_t i = 0; i < B->size2; i++) { cPar.beta_remle_null.push_back(gsl_matrix_get(B, 0, i)); cPar.se_beta_remle_null.push_back(gsl_matrix_get(se_B, 0, i)); @@ -2710,11 +2708,11 @@ void GEMMA::BatchRun(PARAM &cPar) { // calculate and output residuals if (cPar.a_mode == 5) { - gsl_vector *Utu_hat = gsl_vector_alloc(Y->size1); - gsl_vector *Ute_hat = gsl_vector_alloc(Y->size1); - gsl_vector *u_hat = gsl_vector_alloc(Y->size1); - gsl_vector *e_hat = gsl_vector_alloc(Y->size1); - gsl_vector *y_hat = gsl_vector_alloc(Y->size1); + 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); + gsl_vector *e_hat = gsl_vector_safe_alloc(Y->size1); + gsl_vector *y_hat = gsl_vector_safe_alloc(Y->size1); // obtain Utu and Ute gsl_vector_memcpy(y_hat, &UtY_col.vector); @@ -2755,18 +2753,18 @@ void GEMMA::BatchRun(PARAM &cPar) { gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); if (!cPar.file_bfile.empty()) { + // PLINK analysis if (cPar.file_gxe.empty()) { cLmm.AnalyzePlink(U, eval, UtW, &UtY_col.vector, W, - &Y_col.vector); - } else { + &Y_col.vector, cPar.setGWASnps); + } + else { cLmm.AnalyzePlinkGXE(U, eval, UtW, &UtY_col.vector, W, &Y_col.vector, env); } } - // WJA added - else if (!cPar.file_oxford.empty()) { - cLmm.Analyzebgen(U, eval, UtW, &UtY_col.vector, W, &Y_col.vector); - } else { + else { + // BIMBAM analysis if (cPar.file_gxe.empty()) { cLmm.AnalyzeBimbam(U, eval, UtW, &UtY_col.vector, W, &Y_col.vector, cPar.setGWASnps); @@ -2788,8 +2786,6 @@ void GEMMA::BatchRun(PARAM &cPar) { } else { cMvlmm.AnalyzePlinkGXE(U, eval, UtW, UtY, env); } - } else if (!cPar.file_oxford.empty()) { - cMvlmm.Analyzebgen(U, eval, UtW, UtY); } else { if (cPar.file_gxe.empty()) { cMvlmm.AnalyzeBimbam(U, eval, UtW, UtY); @@ -2819,10 +2815,10 @@ void GEMMA::BatchRun(PARAM &cPar) { // BSLMM if (cPar.a_mode == 11 || cPar.a_mode == 12 || cPar.a_mode == 13) { - gsl_vector *y = gsl_vector_alloc(cPar.ni_test); - gsl_matrix *W = gsl_matrix_alloc(y->size, cPar.n_cvt); - gsl_matrix *G = gsl_matrix_alloc(y->size, y->size); - gsl_matrix *UtX = gsl_matrix_alloc(y->size, cPar.ns_test); + gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test); + gsl_matrix *W = gsl_matrix_safe_alloc(y->size, cPar.n_cvt); + gsl_matrix *G = gsl_matrix_safe_alloc(y->size, y->size); + gsl_matrix *UtX = gsl_matrix_safe_alloc(y->size, cPar.ns_test); // set covariates matrix W and phenotype vector y // an intercept should be included in W, @@ -2845,10 +2841,10 @@ void GEMMA::BatchRun(PARAM &cPar) { cBslmm.CopyToParam(cPar); // else, if rho!=1 } else { - gsl_matrix *U = gsl_matrix_alloc(y->size, y->size); - gsl_vector *eval = gsl_vector_alloc(y->size); - gsl_matrix *UtW = gsl_matrix_alloc(y->size, W->size2); - gsl_vector *Uty = gsl_vector_alloc(y->size); + gsl_matrix *U = gsl_matrix_safe_alloc(y->size, y->size); + gsl_vector *eval = gsl_vector_safe_alloc(y->size); + gsl_matrix *UtW = gsl_matrix_safe_alloc(y->size, W->size2); + gsl_vector *Uty = gsl_vector_safe_alloc(y->size); // read relatedness matrix G if (!(cPar.file_kin).empty()) { @@ -2864,7 +2860,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); } @@ -2929,10 +2925,10 @@ void GEMMA::BatchRun(PARAM &cPar) { // BSLMM-DAP if (cPar.a_mode == 14 || cPar.a_mode == 15 || cPar.a_mode == 16) { if (cPar.a_mode == 14) { - gsl_vector *y = gsl_vector_alloc(cPar.ni_test); - gsl_matrix *W = gsl_matrix_alloc(y->size, cPar.n_cvt); - gsl_matrix *G = gsl_matrix_alloc(y->size, y->size); - gsl_matrix *UtX = gsl_matrix_alloc(y->size, cPar.ns_test); + gsl_vector *y = gsl_vector_safe_alloc(cPar.ni_test); + gsl_matrix *W = gsl_matrix_safe_alloc(y->size, cPar.n_cvt); + gsl_matrix *G = gsl_matrix_safe_alloc(y->size, y->size); + gsl_matrix *UtX = gsl_matrix_safe_alloc(y->size, cPar.ns_test); // set covariates matrix W and phenotype vector y // an intercept should be included in W, @@ -2956,10 +2952,10 @@ void GEMMA::BatchRun(PARAM &cPar) { cBslmm.CopyToParam(cPar); // else, if rho!=1 } else { - gsl_matrix *U = gsl_matrix_alloc(y->size, y->size); - gsl_vector *eval = gsl_vector_alloc(y->size); - gsl_matrix *UtW = gsl_matrix_alloc(y->size, W->size2); - gsl_vector *Uty = gsl_vector_alloc(y->size); + gsl_matrix *U = gsl_matrix_safe_alloc(y->size, y->size); + gsl_vector *eval = gsl_vector_safe_alloc(y->size); + gsl_matrix *UtW = gsl_matrix_safe_alloc(y->size, W->size2); + gsl_vector *Uty = gsl_vector_safe_alloc(y->size); // read relatedness matrix G if (!(cPar.file_kin).empty()) { @@ -2975,7 +2971,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); @@ -3090,6 +3086,11 @@ void GEMMA::BatchRun(PARAM &cPar) { return; } +#include "Eigen/Dense" +#if defined(OPENBLAS) && !defined(OPENBLAS_LEGACY) +#include <openblas_config.h> +#endif + void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) { string file_str; file_str = cPar.path_out + "/" + cPar.file_out; @@ -3102,9 +3103,21 @@ void GEMMA::WriteLog(int argc, char **argv, PARAM &cPar) { } outfile << "##" << endl; - outfile << "## GEMMA Version = " << version << endl; - outfile << "## GSL Version = " << GSL_VERSION << endl; - outfile << "## Eigen Version = " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << endl; + outfile << "## GEMMA Version = " << version << endl; + outfile << "## GSL Version = " << GSL_VERSION << endl; + outfile << "## Eigen Version = " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << endl; +#ifdef OPENBLAS + + #ifndef OPENBLAS_LEGACY + outfile << "## OpenBlas =" << OPENBLAS_VERSION << " - " << openblas_get_config() << endl; + outfile << "## arch = " << openblas_get_corename() << endl; + outfile << "## threads = " << openblas_get_num_threads() << endl; + #else + outfile << "## OpenBlas = " << openblas_get_config() << endl; + #endif + string* pStr = new string[4] { "sequential", "threaded", "openmp" }; + outfile << "## parallel type = " << pStr[openblas_get_parallel()] << endl; +#endif outfile << "##" << endl; outfile << "## Command Line Input = "; |