aboutsummaryrefslogtreecommitdiff
path: root/src/gemma.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/gemma.cpp')
-rw-r--r--src/gemma.cpp303
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 = ";
generated by cgit v1.2.3 (git 2.44.0) at 2025-01-19 09:11:06 +0000