diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/debug.cpp | 18 | ||||
| -rw-r--r-- | src/debug.h | 1 | ||||
| -rw-r--r-- | src/fastblas.cpp | 2 | ||||
| -rw-r--r-- | src/gemma.cpp | 193 | ||||
| -rw-r--r-- | src/io.cpp | 92 | ||||
| -rw-r--r-- | src/lmm.cpp | 188 | ||||
| -rw-r--r-- | src/logistic.cpp | 29 | ||||
| -rw-r--r-- | src/mathfunc.cpp | 32 |
8 files changed, 286 insertions, 269 deletions
diff --git a/src/debug.cpp b/src/debug.cpp index eb5c041..4db3daf 100644 --- a/src/debug.cpp +++ b/src/debug.cpp @@ -19,7 +19,7 @@ #include "mathfunc.h" static bool debug_mode = false; -static bool debug_check = false; // check data/algorithms +static bool debug_check = true; // check data/algorithms static bool debug_strict = false; // fail on error static bool debug_quiet = false; static uint debug_issue = 0; // github issues @@ -56,6 +56,22 @@ gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols) { return m; } +/* + Helper function to make sure gsl allocations do their job because + gsl_vector_alloc does not initiatize values (behaviour that changed + in GSL2) we introduced a 'strict mode' by initializing the buffer + with NaNs. This happens when NO-CHECKS is not set + (i.e. -no-checks option). +*/ +gsl_vector *gsl_vector_safe_alloc(size_t n) { + gsl_vector *v = gsl_vector_alloc(n); + enforce_msg(v,"Not enough memory"); // just to be sure when there is no error handler set + if (is_check_mode()) { + gsl_vector_set_all(v, nan("")); + } + return v; +} + // Helper function called by macro validate_K(K, check) void do_validate_K(const gsl_matrix *K, const char *__file, int __line) { if (is_check_mode()) { diff --git a/src/debug.h b/src/debug.h index 7d2197d..c127630 100644 --- a/src/debug.h +++ b/src/debug.h @@ -26,6 +26,7 @@ bool is_issue(uint issue); bool is_legacy_mode(); gsl_matrix *gsl_matrix_safe_alloc(size_t rows,size_t cols); +gsl_vector *gsl_vector_safe_alloc(size_t n); // Validation routines void do_validate_K(const gsl_matrix *K, const char *__file, int __line); diff --git a/src/fastblas.cpp b/src/fastblas.cpp index 35d19c1..de7e27b 100644 --- a/src/fastblas.cpp +++ b/src/fastblas.cpp @@ -197,7 +197,7 @@ void fast_dgemm(const char *TransA, const char *TransB, const double alpha, } else { fast_cblas_dgemm(TransA,TransB,alpha,A,B,beta,C); - #ifndef NDEBUG + #ifdef DISABLE if (is_check_mode()) { // ---- validate with original implementation gsl_matrix *C1 = gsl_matrix_alloc(C->size1,C->size2); diff --git a/src/gemma.cpp b/src/gemma.cpp index 75b1c33..f9a2fc9 100644 --- a/src/gemma.cpp +++ b/src/gemma.cpp @@ -1636,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); @@ -1648,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; @@ -1707,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 @@ -1761,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); @@ -1792,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, @@ -1873,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(); @@ -1918,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); @@ -1928,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); @@ -1971,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); @@ -2029,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, @@ -2082,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); @@ -2218,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 = @@ -2286,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, @@ -2365,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); @@ -2392,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; @@ -2523,19 +2523,19 @@ 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); + 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 @@ -2667,26 +2667,22 @@ void GEMMA::BatchRun(PARAM &cPar) { 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())); @@ -2698,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)); @@ -2709,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); @@ -2816,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, @@ -2842,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()) { @@ -2926,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, @@ -2953,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()) { diff --git a/src/io.cpp b/src/io.cpp index d1c4bc2..8abdeec 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -606,12 +606,12 @@ bool ReadFile_geno(const string &file_geno, const set<string> &setSnps, return false; } - gsl_vector *genotype = gsl_vector_alloc(W->size1); - gsl_vector *genotype_miss = gsl_vector_alloc(W->size1); - gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2); - gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2); - gsl_vector *Wtx = gsl_vector_alloc(W->size2); - gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2); + gsl_vector *genotype = gsl_vector_safe_alloc(W->size1); + gsl_vector *genotype_miss = gsl_vector_safe_alloc(W->size1); + gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2); + gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2); gsl_permutation *pmt = gsl_permutation_alloc(W->size2); gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW); @@ -817,12 +817,12 @@ bool ReadFile_bed(const string &file_bed, const set<string> &setSnps, return false; } - gsl_vector *genotype = gsl_vector_alloc(W->size1); - gsl_vector *genotype_miss = gsl_vector_alloc(W->size1); - gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2); - gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2); - gsl_vector *Wtx = gsl_vector_alloc(W->size2); - gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2); + gsl_vector *genotype = gsl_vector_safe_alloc(W->size1); + gsl_vector *genotype_miss = gsl_vector_safe_alloc(W->size1); + gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2); + gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2); gsl_permutation *pmt = gsl_permutation_alloc(W->size2); gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW); @@ -1366,12 +1366,12 @@ bool BimbamKin(const string file_geno, const set<string> ksnps, bool process_ksnps = ksnps.size(); size_t ni_total = matrix_kin->size1; - gsl_vector *geno = gsl_vector_alloc(ni_total); - gsl_vector *geno_miss = gsl_vector_alloc(ni_total); + gsl_vector *geno = gsl_vector_safe_alloc(ni_total); + gsl_vector *geno_miss = gsl_vector_safe_alloc(ni_total); // Xlarge contains inds x markers const size_t msize = K_BATCH_SIZE; - gsl_matrix *Xlarge = gsl_matrix_alloc(ni_total, msize); + gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_total, msize); enforce_msg(Xlarge, "allocate Xlarge"); gsl_matrix_set_zero(Xlarge); @@ -1506,14 +1506,14 @@ bool PlinkKin(const string &file_bed, vector<int> &indicator_snp, double d, geno_mean, geno_var; size_t ni_total = matrix_kin->size1; - gsl_vector *geno = gsl_vector_alloc(ni_total); + gsl_vector *geno = gsl_vector_safe_alloc(ni_total); size_t ns_test = 0; int n_bit; // Create a large matrix. const size_t msize = K_BATCH_SIZE; - gsl_matrix *Xlarge = gsl_matrix_alloc(ni_total, msize); + gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_total, msize); gsl_matrix_set_zero(Xlarge); // Calculate n_bit and c, the number of bit for each snp. @@ -1649,8 +1649,8 @@ bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv, gsl_matrix_set_zero(K); } - gsl_vector *genotype = gsl_vector_alloc(UtX->size1); - gsl_vector *genotype_miss = gsl_vector_alloc(UtX->size1); + gsl_vector *genotype = gsl_vector_safe_alloc(UtX->size1); + gsl_vector *genotype_miss = gsl_vector_safe_alloc(UtX->size1); double geno, geno_mean; size_t n_miss; @@ -1760,8 +1760,8 @@ bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv, gsl_matrix_set_zero(K); } - gsl_vector *genotype = gsl_vector_alloc(ni_test); - gsl_vector *genotype_miss = gsl_vector_alloc(ni_test); + gsl_vector *genotype = gsl_vector_safe_alloc(ni_test); + gsl_vector *genotype_miss = gsl_vector_safe_alloc(ni_test); double geno, geno_mean; size_t n_miss; @@ -1879,7 +1879,7 @@ bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv, gsl_matrix_set_zero(K); } - gsl_vector *genotype = gsl_vector_alloc(UtX->size1); + gsl_vector *genotype = gsl_vector_safe_alloc(UtX->size1); double geno, geno_mean; size_t n_miss; @@ -2015,7 +2015,7 @@ bool ReadFile_bed(const string &file_bed, vector<int> &indicator_idv, gsl_matrix_set_zero(K); } - gsl_vector *genotype = gsl_vector_alloc(ni_test); + gsl_vector *genotype = gsl_vector_safe_alloc(ni_test); double geno, geno_mean; size_t n_miss; @@ -2658,13 +2658,13 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps, double d, geno_mean, geno_var; size_t ni_test = matrix_kin->size1; - gsl_vector *geno = gsl_vector_alloc(ni_test); - gsl_vector *geno_miss = gsl_vector_alloc(ni_test); + gsl_vector *geno = gsl_vector_safe_alloc(ni_test); + gsl_vector *geno_miss = gsl_vector_safe_alloc(ni_test); - gsl_vector *Wtx = gsl_vector_alloc(W->size2); - gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2); - gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2); - gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2); + gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2); + gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2); gsl_permutation *pmt = gsl_permutation_alloc(W->size2); gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW); @@ -2681,7 +2681,7 @@ bool BimbamKinUncentered(const string &file_geno, const set<string> ksnps, // Create a large matrix. const size_t msize = K_BATCH_SIZE; - gsl_matrix *Xlarge = gsl_matrix_alloc(ni_test, msize * n_vc); + gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_test, msize * n_vc); gsl_matrix_set_zero(Xlarge); size_t ns_test = 0; @@ -2850,12 +2850,12 @@ bool PlinkKin(const string &file_bed, const int display_pace, size_t ni_test = matrix_kin->size1; size_t ni_total = indicator_idv.size(); - gsl_vector *geno = gsl_vector_alloc(ni_test); + gsl_vector *geno = gsl_vector_safe_alloc(ni_test); - gsl_vector *Wtx = gsl_vector_alloc(W->size2); - gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2); - gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2); - gsl_vector *WtWiWtx = gsl_vector_alloc(W->size2); + gsl_vector *Wtx = gsl_vector_safe_alloc(W->size2); + gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_vector *WtWiWtx = gsl_vector_safe_alloc(W->size2); gsl_permutation *pmt = gsl_permutation_alloc(W->size2); gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW); @@ -2875,7 +2875,7 @@ bool PlinkKin(const string &file_bed, const int display_pace, // Create a large matrix. const size_t msize = K_BATCH_SIZE; - gsl_matrix *Xlarge = gsl_matrix_alloc(ni_test, msize * n_vc); + gsl_matrix *Xlarge = gsl_matrix_safe_alloc(ni_test, msize * n_vc); gsl_matrix_set_zero(Xlarge); // Calculate n_bit and c, the number of bit for each SNP. @@ -3074,8 +3074,8 @@ bool MFILEKin(const size_t mfile_mode, const string &file_mfile, string file_name; - gsl_matrix *kin_tmp = gsl_matrix_alloc(matrix_kin->size1, matrix_kin->size2); - gsl_vector *ns_tmp = gsl_vector_alloc(vector_ns->size); + gsl_matrix *kin_tmp = gsl_matrix_safe_alloc(matrix_kin->size1, matrix_kin->size2); + gsl_vector *ns_tmp = gsl_vector_safe_alloc(vector_ns->size); size_t l = 0; double d; @@ -3843,7 +3843,7 @@ void ReadFile_study(const string &file_study, gsl_matrix *Vq_mat, string sfile = file_study + ".size.txt"; string qfile = file_study + ".q.txt"; - gsl_vector *s = gsl_vector_alloc(s_vec->size + 1); + gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1); ReadFile_matrix(Vqfile, Vq_mat); ReadFile_vector(sfile, s); @@ -3868,7 +3868,7 @@ void ReadFile_ref(const string &file_ref, gsl_matrix *S_mat, string sfile = file_ref + ".size.txt"; string Sfile = file_ref + ".S.txt"; - gsl_vector *s = gsl_vector_alloc(s_vec->size + 1); + gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1); ReadFile_vector(sfile, s); ReadFile_matrix(Sfile, S_mat, Svar_mat); @@ -3894,9 +3894,9 @@ void ReadFile_mstudy(const string &file_mstudy, gsl_matrix *Vq_mat, gsl_vector_set_zero(s_vec); ni = 0; - gsl_matrix *Vq_sub = gsl_matrix_alloc(Vq_mat->size1, Vq_mat->size2); - gsl_vector *q_sub = gsl_vector_alloc(q_vec->size); - gsl_vector *s = gsl_vector_alloc(s_vec->size + 1); + gsl_matrix *Vq_sub = gsl_matrix_safe_alloc(Vq_mat->size1, Vq_mat->size2); + gsl_vector *q_sub = gsl_vector_safe_alloc(q_vec->size); + gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1); igzstream infile(file_mstudy.c_str(), igzstream::in); if (!infile) { @@ -3985,9 +3985,9 @@ void ReadFile_mref(const string &file_mref, gsl_matrix *S_mat, gsl_vector_set_zero(s_vec); ni = 0; - gsl_matrix *S_sub = gsl_matrix_alloc(S_mat->size1, S_mat->size2); - gsl_matrix *Svar_sub = gsl_matrix_alloc(Svar_mat->size1, Svar_mat->size2); - gsl_vector *s = gsl_vector_alloc(s_vec->size + 1); + gsl_matrix *S_sub = gsl_matrix_safe_alloc(S_mat->size1, S_mat->size2); + gsl_matrix *Svar_sub = gsl_matrix_safe_alloc(Svar_mat->size1, Svar_mat->size2); + gsl_vector *s = gsl_vector_safe_alloc(s_vec->size + 1); igzstream infile(file_mref.c_str(), igzstream::in); if (!infile) { diff --git a/src/lmm.cpp b/src/lmm.cpp index 7b32330..e80cd76 100644 --- a/src/lmm.cpp +++ b/src/lmm.cpp @@ -363,9 +363,9 @@ double LogL_f(double l, void *params) { double f = 0.0, logdet_h = 0.0, d; size_t index_yy; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -413,11 +413,11 @@ double LogL_dev1(double l, void *params) { double dev1 = 0.0, trace_Hi = 0.0; size_t index_yy; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -476,13 +476,13 @@ double LogL_dev2(double l, void *params) { double dev2 = 0.0, trace_Hi = 0.0, trace_HiHi = 0.0; size_t index_yy; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -553,13 +553,13 @@ void LogL_dev12(double l, void *params, double *dev1, double *dev2) { double trace_Hi = 0.0, trace_HiHi = 0.0; size_t index_yy; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -636,10 +636,10 @@ double LogRL_f(double l, void *params) { double f = 0.0, logdet_h = 0.0, logdet_hiw = 0.0, d; size_t index_ww; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *Iab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *Iab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -701,11 +701,11 @@ double LogRL_dev1(double l, void *params) { double dev1 = 0.0, trace_Hi = 0.0; size_t index_ww; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -777,13 +777,13 @@ double LogRL_dev2(double l, void *params) { double dev2 = 0.0, trace_Hi = 0.0, trace_HiHi = 0.0; size_t index_ww; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -867,13 +867,13 @@ void LogRL_dev12(double l, void *params, double *dev1, double *dev2) { double trace_Hi = 0.0, trace_HiHi = 0.0; size_t index_ww; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_matrix *PPPab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *HiHiHi_eval = gsl_vector_alloc((p->eval)->size); - gsl_vector *v_temp = gsl_vector_alloc((p->eval)->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_matrix *PPPab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *HiHiHi_eval = gsl_vector_safe_alloc((p->eval)->size); + gsl_vector *v_temp = gsl_vector_safe_alloc((p->eval)->size); gsl_vector_memcpy(v_temp, p->eval); gsl_vector_scale(v_temp, l); @@ -947,9 +947,9 @@ void LMM::CalcRLWald(const double &l, const FUNC_PARAM ¶ms, double &beta, int df = (int)ni_test - (int)n_cvt - 1; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc(params.eval->size); - gsl_vector *v_temp = gsl_vector_alloc(params.eval->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc(params.eval->size); + gsl_vector *v_temp = gsl_vector_safe_alloc(params.eval->size); gsl_vector_memcpy(v_temp, params.eval); gsl_vector_scale(v_temp, l); @@ -989,9 +989,9 @@ void LMM::CalcRLScore(const double &l, const FUNC_PARAM ¶ms, double &beta, int df = (int)ni_test - (int)n_cvt - 1; - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc(params.eval->size); - gsl_vector *v_temp = gsl_vector_alloc(params.eval->size); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc(params.eval->size); + gsl_vector *v_temp = gsl_vector_safe_alloc(params.eval->size); gsl_vector_memcpy(v_temp, params.eval); gsl_vector_scale(v_temp, l); @@ -1030,7 +1030,7 @@ void CalcUab(const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) { size_t index_ab; size_t n_cvt = UtW->size2; - gsl_vector *u_a = gsl_vector_alloc(Uty->size); + gsl_vector *u_a = gsl_vector_safe_alloc(Uty->size); for (size_t a = 1; a <= n_cvt + 2; ++a) { if (a == n_cvt + 1) { @@ -1096,8 +1096,8 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) { size_t n_cvt = W->size2; double d; - gsl_vector *v_a = gsl_vector_alloc(y->size); - gsl_vector *v_b = gsl_vector_alloc(y->size); + gsl_vector *v_a = gsl_vector_safe_alloc(y->size); + gsl_vector *v_b = gsl_vector_safe_alloc(y->size); for (size_t a = 1; a <= n_cvt + 2; ++a) { if (a == n_cvt + 1) { @@ -1141,7 +1141,7 @@ void Calcab(const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x, size_t n_cvt = W->size2; double d; - gsl_vector *v_b = gsl_vector_alloc(y->size); + gsl_vector *v_b = gsl_vector_safe_alloc(y->size); for (size_t b = 1; b <= n_cvt + 2; ++b) { index_ab = GetabIndex(n_cvt + 1, b, n_cvt); @@ -1188,10 +1188,10 @@ void LMM::AnalyzeGene(const gsl_matrix *U, const gsl_vector *eval, // Calculate basic quantities. size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2; - gsl_vector *y = gsl_vector_alloc(U->size1); - gsl_vector *Uty = gsl_vector_alloc(U->size2); - gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index); - gsl_vector *ab = gsl_vector_alloc(n_index); + gsl_vector *y = gsl_vector_safe_alloc(U->size1); + gsl_vector *Uty = gsl_vector_safe_alloc(U->size2); + gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index); + gsl_vector *ab = gsl_vector_safe_alloc(n_index); // Header. getline(infile, line); @@ -1289,17 +1289,17 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp, const size_t inds = U->size1; enforce(inds == ni_test); - gsl_vector *x = gsl_vector_alloc(inds); // #inds - gsl_vector *x_miss = gsl_vector_alloc(inds); - gsl_vector *Utx = gsl_vector_alloc(U->size2); - gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index); - gsl_vector *ab = gsl_vector_alloc(n_index); + gsl_vector *x = gsl_vector_safe_alloc(inds); // #inds + gsl_vector *x_miss = gsl_vector_safe_alloc(inds); + gsl_vector *Utx = gsl_vector_safe_alloc(U->size2); + gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index); + gsl_vector *ab = gsl_vector_safe_alloc(n_index); // Create a large matrix with LMM_BATCH_SIZE columns for batched processing // const size_t msize=(process_gwasnps ? 1 : LMM_BATCH_SIZE); const size_t msize = LMM_BATCH_SIZE; - gsl_matrix *Xlarge = gsl_matrix_alloc(inds, msize); - gsl_matrix *UtXlarge = gsl_matrix_alloc(inds, msize); + gsl_matrix *Xlarge = gsl_matrix_safe_alloc(inds, msize); + gsl_matrix *UtXlarge = gsl_matrix_safe_alloc(inds, msize); enforce_msg(Xlarge && UtXlarge, "Xlarge memory check"); // just to be sure enforce(Xlarge->size1 == inds); gsl_matrix_set_zero(Xlarge); @@ -1573,10 +1573,10 @@ void MatrixCalcLR(const gsl_matrix *U, const gsl_matrix *UtX, vector<pair<size_t, double>> &pos_loglr) { double logl_H0, logl_H1, log_lr, lambda0, lambda1; - gsl_vector *w = gsl_vector_alloc(Uty->size); - gsl_matrix *Utw = gsl_matrix_alloc(Uty->size, 1); - gsl_matrix *Uab = gsl_matrix_alloc(Uty->size, 6); - gsl_vector *ab = gsl_vector_alloc(6); + gsl_vector *w = gsl_vector_safe_alloc(Uty->size); + gsl_matrix *Utw = gsl_matrix_safe_alloc(Uty->size, 1); + gsl_matrix *Uab = gsl_matrix_safe_alloc(Uty->size, 6); + gsl_vector *ab = gsl_vector_safe_alloc(6); gsl_vector_set_zero(ab); gsl_vector_set_all(w, 1.0); @@ -1781,8 +1781,8 @@ void CalcLambda(const char func_name, const gsl_vector *eval, size_t n_cvt = UtW->size2, ni_test = UtW->size1; size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2; - gsl_matrix *Uab = gsl_matrix_alloc(ni_test, n_index); - gsl_vector *ab = gsl_vector_alloc(n_index); + gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index); + gsl_vector *ab = gsl_vector_safe_alloc(n_index); gsl_matrix_set_zero(Uab); CalcUab(UtW, Uty, Uab); @@ -1804,8 +1804,8 @@ void CalcPve(const gsl_vector *eval, const gsl_matrix *UtW, size_t n_cvt = UtW->size2, ni_test = UtW->size1; size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2; - gsl_matrix *Uab = gsl_matrix_alloc(ni_test, n_index); - gsl_vector *ab = gsl_vector_alloc(n_index); + gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index); + gsl_vector *ab = gsl_vector_safe_alloc(n_index); gsl_matrix_set_zero(Uab); CalcUab(UtW, Uty, Uab); @@ -1831,15 +1831,15 @@ void CalcLmmVgVeBeta(const gsl_vector *eval, const gsl_matrix *UtW, size_t n_cvt = UtW->size2, ni_test = UtW->size1; size_t n_index = (n_cvt + 2 + 1) * (n_cvt + 2) / 2; - gsl_matrix *Uab = gsl_matrix_alloc(ni_test, n_index); - gsl_vector *ab = gsl_vector_alloc(n_index); - gsl_matrix *Pab = gsl_matrix_alloc(n_cvt + 2, n_index); - gsl_vector *Hi_eval = gsl_vector_alloc(eval->size); - gsl_vector *v_temp = gsl_vector_alloc(eval->size); - gsl_matrix *HiW = gsl_matrix_alloc(eval->size, UtW->size2); - gsl_matrix *WHiW = gsl_matrix_alloc(UtW->size2, UtW->size2); - gsl_vector *WHiy = gsl_vector_alloc(UtW->size2); - gsl_matrix *Vbeta = gsl_matrix_alloc(UtW->size2, UtW->size2); + gsl_matrix *Uab = gsl_matrix_safe_alloc(ni_test, n_index); + gsl_vector *ab = gsl_vector_safe_alloc(n_index); + gsl_matrix *Pab = gsl_matrix_safe_alloc(n_cvt + 2, n_index); + gsl_vector *Hi_eval = gsl_vector_safe_alloc(eval->size); + gsl_vector *v_temp = gsl_vector_safe_alloc(eval->size); + gsl_matrix *HiW = gsl_matrix_safe_alloc(eval->size, UtW->size2); + gsl_matrix *WHiW = gsl_matrix_safe_alloc(UtW->size2, UtW->size2); + gsl_vector *WHiy = gsl_vector_safe_alloc(UtW->size2); + gsl_matrix *Vbeta = gsl_matrix_safe_alloc(UtW->size2, UtW->size2); gsl_matrix_set_zero(Uab); CalcUab(UtW, Uty, Uab); @@ -1921,13 +1921,13 @@ void LMM::AnalyzeBimbamGXE(const gsl_matrix *U, const gsl_vector *eval, // Calculate basic quantities. size_t n_index = (n_cvt + 2 + 2 + 1) * (n_cvt + 2 + 2) / 2; - gsl_vector *x = gsl_vector_alloc(U->size1); - gsl_vector *x_miss = gsl_vector_alloc(U->size1); - gsl_vector *Utx = gsl_vector_alloc(U->size2); - gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index); - gsl_vector *ab = gsl_vector_alloc(n_index); + gsl_vector *x = gsl_vector_safe_alloc(U->size1); + gsl_vector *x_miss = gsl_vector_safe_alloc(U->size1); + gsl_vector *Utx = gsl_vector_safe_alloc(U->size2); + gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index); + gsl_vector *ab = gsl_vector_safe_alloc(n_index); - gsl_matrix *UtW_expand = gsl_matrix_alloc(U->size1, UtW->size2 + 2); + gsl_matrix *UtW_expand = gsl_matrix_safe_alloc(U->size1, UtW->size2 + 2); gsl_matrix_view UtW_expand_mat = gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2); gsl_matrix_memcpy(&UtW_expand_mat.matrix, UtW); @@ -2071,12 +2071,12 @@ void LMM::AnalyzePlinkGXE(const gsl_matrix *U, const gsl_vector *eval, // Calculate basic quantities. size_t n_index = (n_cvt + 2 + 2 + 1) * (n_cvt + 2 + 2) / 2; - gsl_vector *x = gsl_vector_alloc(U->size1); - gsl_vector *Utx = gsl_vector_alloc(U->size2); - gsl_matrix *Uab = gsl_matrix_alloc(U->size2, n_index); - gsl_vector *ab = gsl_vector_alloc(n_index); + gsl_vector *x = gsl_vector_safe_alloc(U->size1); + gsl_vector *Utx = gsl_vector_safe_alloc(U->size2); + gsl_matrix *Uab = gsl_matrix_safe_alloc(U->size2, n_index); + gsl_vector *ab = gsl_vector_safe_alloc(n_index); - gsl_matrix *UtW_expand = gsl_matrix_alloc(U->size1, UtW->size2 + 2); + gsl_matrix *UtW_expand = gsl_matrix_safe_alloc(U->size1, UtW->size2 + 2); gsl_matrix_view UtW_expand_mat = gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2); gsl_matrix_memcpy(&UtW_expand_mat.matrix, UtW); diff --git a/src/logistic.cpp b/src/logistic.cpp index 2dd0402..a936682 100644 --- a/src/logistic.cpp +++ b/src/logistic.cpp @@ -7,6 +7,7 @@ #include <stdio.h> #include "logistic.h" +#include "debug.h" // I need to bundle all the data that goes to the function to optimze // together. @@ -135,7 +136,7 @@ void wgsl_mixed_optim_hessian(const gsl_vector *beta, void *params, int K = p->X->size2; int Kc = p->Xc->size2; int npar = beta->size; - gsl_vector *gn = gsl_vector_alloc(npar); // gn + gsl_vector *gn = gsl_vector_safe_alloc(npar); // gn // Intitialize Hessian out necessary ??? gsl_matrix_set_zero(out); @@ -226,11 +227,11 @@ int logistic_mixed_fit(gsl_vector *beta, gsl_matrix_int *X, // Initial fit. mLogLik = wgsl_mixed_optim_f(beta, &p); - gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix. - gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. + gsl_matrix *myH = gsl_matrix_safe_alloc(npar, npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_safe_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. - gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. + gsl_vector *myG = gsl_vector_safe_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_safe_alloc(npar); // tau for QR. for (iter = 0; iter < 100; iter++) { wgsl_mixed_optim_hessian(beta, &p, myH); // Calculate Hessian. @@ -456,11 +457,11 @@ int logistic_cat_fit(gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev, // Initial fit. mLogLik = wgsl_cat_optim_f(beta, &p); - gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix. - gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. + gsl_matrix *myH = gsl_matrix_safe_alloc(npar, npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_safe_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. - gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. + gsl_vector *myG = gsl_vector_safe_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_safe_alloc(npar); // tau for QR. for (iter = 0; iter < 100; iter++) { wgsl_cat_optim_hessian(beta, &p, myH); // Calculate Hessian. @@ -596,7 +597,7 @@ void wgsl_cont_optim_hessian(const gsl_vector *beta, void *params, int n = p->y->size; int Kc = p->Xc->size2; int npar = beta->size; - gsl_vector *gn = gsl_vector_alloc(npar); // gn. + gsl_vector *gn = gsl_vector_safe_alloc(npar); // gn. // Intitialize Hessian out necessary ??? @@ -673,11 +674,11 @@ int logistic_cont_fit(gsl_vector *beta, // Initial fit. mLogLik = wgsl_cont_optim_f(beta, &p); - gsl_matrix *myH = gsl_matrix_alloc(npar, npar); // Hessian matrix. - gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. + gsl_matrix *myH = gsl_matrix_safe_alloc(npar, npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_safe_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. - gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. + gsl_vector *myG = gsl_vector_safe_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_safe_alloc(npar); // tau for QR. for (iter = 0; iter < 100; iter++) { wgsl_cont_optim_hessian(beta, &p, myH); // Calculate Hessian. diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp index f55f3c6..e7dff73 100644 --- a/src/mathfunc.cpp +++ b/src/mathfunc.cpp @@ -80,8 +80,8 @@ double VectorVar(const gsl_vector *v) { // Center the matrix G. void CenterMatrix(gsl_matrix *G) { double d; - gsl_vector *w = gsl_vector_alloc(G->size1); - gsl_vector *Gw = gsl_vector_alloc(G->size1); + gsl_vector *w = gsl_vector_safe_alloc(G->size1); + gsl_vector *Gw = gsl_vector_safe_alloc(G->size1); gsl_vector_set_all(w, 1.0); gsl_blas_dgemv(CblasNoTrans, 1.0, G, w, 0.0, Gw); @@ -105,7 +105,7 @@ void CenterMatrix(gsl_matrix *G) { // Center the matrix G. void CenterMatrix(gsl_matrix *G, const gsl_vector *w) { double d, wtw; - gsl_vector *Gw = gsl_vector_alloc(G->size1); + gsl_vector *Gw = gsl_vector_safe_alloc(G->size1); gsl_blas_ddot(w, w, &wtw); gsl_blas_dgemv(CblasNoTrans, 1.0, G, w, 0.0, Gw); @@ -127,12 +127,12 @@ void CenterMatrix(gsl_matrix *G, const gsl_vector *w) { // Center the matrix G. void CenterMatrix(gsl_matrix *G, const gsl_matrix *W) { - gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2); - gsl_matrix *WtWi = gsl_matrix_alloc(W->size2, W->size2); - gsl_matrix *WtWiWt = gsl_matrix_alloc(W->size2, G->size1); - gsl_matrix *GW = gsl_matrix_alloc(G->size1, W->size2); - gsl_matrix *WtGW = gsl_matrix_alloc(W->size2, W->size2); - gsl_matrix *Gtmp = gsl_matrix_alloc(G->size1, G->size1); + gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_matrix *WtWi = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_matrix *WtWiWt = gsl_matrix_safe_alloc(W->size2, G->size1); + gsl_matrix *GW = gsl_matrix_safe_alloc(G->size1, W->size2); + gsl_matrix *WtGW = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_matrix *Gtmp = gsl_matrix_safe_alloc(G->size1, G->size1); gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW); @@ -227,7 +227,7 @@ bool isMatrixSymmetric(const gsl_matrix *G) { bool isMatrixPositiveDefinite(const gsl_matrix *G) { enforce(G->size1 == G->size2); - auto G2 = gsl_matrix_alloc(G->size1, G->size2); + auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2); enforce_gsl(gsl_matrix_memcpy(G2,G)); auto handler = gsl_set_error_handler_off(); #if GSL_MAJOR_VERSION >= 2 && GSL_MINOR_VERSION >= 3 @@ -242,11 +242,11 @@ bool isMatrixPositiveDefinite(const gsl_matrix *G) { gsl_vector *getEigenValues(const gsl_matrix *G) { enforce(G->size1 == G->size2); - auto G2 = gsl_matrix_alloc(G->size1, G->size2); + auto G2 = gsl_matrix_safe_alloc(G->size1, G->size2); enforce_gsl(gsl_matrix_memcpy(G2,G)); auto eworkspace = gsl_eigen_symm_alloc(G->size1); enforce(eworkspace); - gsl_vector *eigenvalues = gsl_vector_alloc(G->size1); + gsl_vector *eigenvalues = gsl_vector_safe_alloc(G->size1); enforce_gsl(gsl_eigen_symm(G2, eigenvalues, eworkspace)); gsl_eigen_symm_free(eworkspace); gsl_matrix_free(G2); @@ -345,9 +345,9 @@ double CenterVector(gsl_vector *y) { // Center the vector y. void CenterVector(gsl_vector *y, const gsl_matrix *W) { - gsl_matrix *WtW = gsl_matrix_alloc(W->size2, W->size2); - gsl_vector *Wty = gsl_vector_alloc(W->size2); - gsl_vector *WtWiWty = gsl_vector_alloc(W->size2); + gsl_matrix *WtW = gsl_matrix_safe_alloc(W->size2, W->size2); + gsl_vector *Wty = gsl_vector_safe_alloc(W->size2); + gsl_vector *WtWiWty = gsl_vector_safe_alloc(W->size2); gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW); gsl_blas_dgemv(CblasTrans, 1.0, W, y, 0.0, Wty); @@ -387,7 +387,7 @@ void StandardizeVector(gsl_vector *y) { // Calculate UtX. void CalcUtX(const gsl_matrix *U, gsl_matrix *UtX) { - gsl_matrix *X = gsl_matrix_alloc(UtX->size1, UtX->size2); + gsl_matrix *X = gsl_matrix_safe_alloc(UtX->size1, UtX->size2); gsl_matrix_memcpy(X, UtX); fast_dgemm("T", "N", 1.0, U, X, 0.0, UtX); gsl_matrix_free(X); |